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aeroelastic  beam  bending  rigidity  about  local  £1  axis  =  JJ  E^dA;  also,  rigid- 

-jV" 

body  mass  principal  moment  of  inertia  for  mass  center  about 
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SUMMARY 


The  Rotorcraft  Dynamics  Division,  Aeroflightdynamics  Directorate,  U.S.  Army  Avia¬ 
tion  Research  and  Technology  Activity  (AVSCOM)  has  developed  the  General  Rotorcraft 
Aeromechanical  Stability  Program  (GRASP)  to  calculate  aeroelastic  stability  for  rotorcraft 
in  hovering  flight,  vertical  flight,  and  ground  contact  conditions.  In  this  report,  GRASP  is 
described  in  terms  of  its  capabilities  and  the  philosophy  behind  its  modeling.  The  equa¬ 
tions  of  motion  that  govern  the  physical  system  are  described,  as  well  as  the  analytical 
approximations  used  to  derive  the  equations.  These  equations  include  the  kinematical 
equation,  the  element  equations,  and  the  constraint  equations.  In  addition,  the  solution 
procedures  used  by  GRASP  are  described.  \ 

GRASP  is  capable  of  treating  the  nonlinear  static  and  linearized  dynamic  behavior 
of  structures  represented  by  arbitrary  collections  of  rigid-body  and  beam  elements.  These 
elements  may  be  connected  in  an  arbitrary  fashion,  and  are  permitted  to  have  large  rel¬ 
ative  motions.  The  main  limitation  of  this  analysis  is  that  periodic  coefficient  effects  are 
not  treated,  restricting  rotorcraft  flight  conditions  to  hover,  axial  flight,  and  ground  con¬ 
tact.  Instead  of  following,  the  methods  employed  in  other  rotorcraft  programs,  GRASP  is 
designed  to  be  a  hybrid  of  the  finite-element  method  and  the  multibody  methods  used  in 
spacecraft  analyses.  GRASP  differs  from  traditional  finite-element  programs  by  allowing 
multiple  levels  of  substructures  in  which  the  substructures  can  move  and/or  rotate  relative 
to  others  with  no  small-angle  approximations.  This  capability  facilitates  the  modeling  of 
rotorcraft  structures,  including  the  rotating/nonrotating  interface  and  the  details  of  the 
blade/root  kinematics  for  various  rotor  types.  GRASP  differs  from  traditional  multibody 
programs  by  considering  aeroelastic  effects,  including  inflow  dynamics  (simple  unsteady 
aerodynamics)  and  nonlinear  aerodynamic  coefficients. 
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1.  INTRODUCTION 


Previcms  helicopter  aeroelastic  stability  programs  have  suffered  from  significant  re¬ 
strictions.  The  General  Rotorcraft  Aeromechanical  Stability  Program  has  been  developed 
using  a  modern  approach  which  overcomes  these  limitations. 


1.1.  Background 

In  early  efforts  made  to  calculate  the  aeroelastic  stability  of  hingeless  helicopter  ro¬ 
tor  blades,  it  was  common  practice  to  make  use  of  simple  physical  models  ( e.g .,  spring- 
restrained,  centrally- hinged,  rigid  blades  (ref.  1)).  Later  work  treated  configurations  that 
were  somewhat  more  complex,  and  included  models  of  elastic  blades  (ref.  2),  body  degrees 
of  freedom,  and  inflow  dynamics  (ref.  3).  These  simple  approaches  to  rotorcraft  aeroelastic 
stability  calculations  have  been  very  valuable  for  gaining  physical  insight  into  many  com¬ 
plicated  phenomena  (e.g.,  coupled  rotor-fuselage  stability).  They  all  are,  however,  based 
on  a  single  physical  model,  and  therefore  are  of  limited  value  when  more  realistic  rotorcraft 
configurations  must  be  analyzed. 

Because  of  the  complex  couplings  inherent  in  a  bent  and  twisted  beam,  the  calculation 
of  aeroelastic  stability  is  particularly  important  in  the  analysis  of  rotor  blades  having 
cantilever  root  boundary  conditions  (e.g.,  hingeless  and  bearingless  rotors).  In  bearingless 
rotors,  the  blade/root  kinematics  demand  a  great  deal  of  modeling  flexibility  because 
individual  blade  designs  tend  to  have  widely  varying  configurations.  The  FLAIR  program 
(refs.  4,  5,  and  6)  is  able  to  perform  this  type  of  aeroelastic  stability  calculation,  but  is 
limited  to  a  configuration  that  has  a  rigid  blade,  a  uniform  flexbeam,  linear  aerodynamics, 
static  induced  velocity,  and  several  different  blade/root  configurations.  While  FLAIR  is 
currently  being  used  in  the  rotorcraft  community,  it  lacks  the  flexibility  and  generality 
necessary  for  it  to  be  considered  general-purpose  analysis. 

For  analysis  of  problems  involving  complete  rotorcraft,  there  exist  large  helicopter  sim¬ 
ulation  programs  such  as  C-81  (ref.  7)  and  G400  (ref.  8).  These  programs  were  designed 
primarily  for  time-history  analysis  of  rotorcraft  behavior  in  forward  flight  rathei  than  for 
aeromechanical  stability.  Despite  their  generality  and  complexity,  these  programs  have 
limitations  (primarily  related  to  aerodynamics)  which  are  pointed  out  by  Johnson  (ref.  9) 
in  his  discussion  of  these  and  other  large  rotorcraft  programs.  While  the  CAMRAD  pro¬ 
gram  overcomes  many  of  these  limitations,  all  of  these  programs  (including  CAMRAD)  are 
restricted  to  a  fixed  number  of  physical  models,  and  lack  the  modeling  flexibility  needed  to 
deal  with  a  wide  variety  of  blade/root  geometries.  Many  of  these  programs  rely  on  results, 
such  as  a  set  of  modes,  from  other  programs.  This  approach  may  present  an  assortment  of 
modeling  difficulties,  especially  for  bcaringless  rotor  blades.  In  particular,  the  mathemati 
cal  and  physical  consistency  of  a  combined  approach  is  seldom  examined,  and  the  physical 
bases  of  the  individual  programs  are  likely  to  not  be  consistent.  Furthermore,  in  stability 
analyses  a  nonlinear  static  equilibrium  solution  is  needed  about  which  to  linearize  —  an 
important  consideration  which  most  of  the  earlier  simulation  programs  do  not  address. 
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Therefore,  it  is  important  that  a  code  be  developed  in  which  blade  structural  dynam¬ 
ics,  isolated  blade  stability,  and  isolated  rotor  stability,  as  well  as  coupled  rotor /airframe 
stability,  can  all  be  treated  under  a  consistent  set  of  physical  assumptions. 

Dynamic  coupling  programs,  such  as  DYSCO  (ref.  10),  which  have  a  high  degree 
of  generality,  allow  coupling  of  discrete  component  models  and/or  modal  representations 
of  flexible  structures.  While  DYSCO  has  a  very  powerful,  executive-driven  system,  it 
currently  cannot  treat  the  aeroelastic  behavior  of  bearinglcss  rotor  systems  undergoing 
geometrically  nonlinear  deformation.  The  problem  is  that  it  lacks  a  sufficiently  general 
element  in  its  element  library. 

Several  recent  implementations  that  apply  the  finite-element  method  to  rotorcraft 
problems  (refs.  11,  12,  and  13)  are  not  abie  to  overcome  these  limitations  because  their 
physical  models  are  limited  to  a  single  configuration.  Simply  breaking  a  rotating  beam 
into  a  number  of  finite  elements  yields  nothing  more  than  a  discretized  rotating  beam. 
This  approach  does  not  meet  the  requirement  that  the  beam  be  coupled  with  an  airframe, 
or  model  blade/root  kinematics  of  an  arbitrary  configuration.  The  classical  finite-element 
method  is  based  on  the  breaking  up  of  a  single  structure  (»’.e.,  a  beam,  plate,  or  shell) 
into  an  arbitrary  number  of  elements  and  expanding  the  appropriate  field  variables  into 
polynomial  shape  functions.  This  approach  by  itself  also  lacks  the  flexibility  to  deal  with 
truly  arbitrary  rotorcraft  configurations  because  a  helicopter  is  a  system  of  structural 
components,  some  of  which  may  be  rotating  and/or  translating  relative  to  one  another. 
Because  of  this,  rotorcraft  are  actually  more  akin  to  the  multibody  systems  (refs.  14  and 
15)  encountered  in  spacecraft  problems.  Unfortunately,  few  multibody  programs  possess 
the  capability  to  deal  with  flexible  components,  and  none  have  the  capability  to  deal  with 
aeroelastic  phenomena  since  they  were  developed  primarily  for  spacecraft  applications. 

All  previous  attempts  at  modeling  rotorcraft  problems  have  incorporated  certain  re¬ 
strictions  that  are  undesirable  in  a  truly  general-purpose  program.  General-purpose  codes 
that  are  currently  under  development,  or  will  be  developed  in  the  future,  should  over¬ 
come  the  major  shortcomings  of  existing  aeroelastic  analyses.  Consider,  for  example,  the 
following  typical  restrictions: 

The  first  is  a  restriction  to  linear,  small  displacement  approximations  of  beam  elastic 
deformation.  This  restriction  is  unacceptable  in  a  general-purpose  rotorcraft  program 
because  the  rotor  blade  aeroelastic  problem,  especially  for  hingeless  and  bearingless  rotor 
blades,  has  been  conclusively  shown  to  be  a  nonlinear  problem.  A  consistent  approach 
based  on  nonlinear  kinematics  is  required  for  these  configurations. 

The  second  is  a  restriction  to  elastic  blade  models  with  ordering  schemes,  second- 
degree  nonlinearity,  or  “moderate”  rotations.  These  approximations  are  undesirable  be¬ 
cause  the  governing  equations  often  have  to  be  augmented  with  certain  higher  order  terms 
if  the  values  of  certain  structural  properties  are  not  within  some  nominal  range  (Rosen 
and  Friedmann  (ref.  16)).  Therefore,  in  a  general  purpose  analysis,  the  higher  order  terms 
must  be  present.  Ordering  schemes,  while  still  a  valuable  tool  when  used  in  special  purpose 
codes  and  codes  where  accuracy  is  a  secondary  consideration,  are  neither  necessary  nor 
desirable  in  a  general-purpose  context.  Furthermore,  a  bearingiess  rotor  fiexbeam  must 
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undergo  deformation-induced  rotations  of  the  order  of  the  collective  pitch  angle  —  a  ro¬ 
tation  too  large  to  be  classified  as  “moderate.”  Thus,  bearingless  rotor  problems  demand 
a  large- deflection  analysis  without  artificial  restrictions  on  rotations  due  to  deformation, 
the  degree  of  nonlinearity,  or  the  values  of  blade  properties. 

The  third  restriction  is  to  a  fixed  number  (usually  one)  of  configurations  ( e.g .,  isolated 
hingeless  blade  or  coupled  bearingless  rotor  and  body  or  a  single  blade/root  configuration). 
This  restriction  is  unacceptable  in  a  general-purpose  code  because  the  intent  of  such  a  code 
is  to  analyze  different  types  of  configurations  with  a  single,  consistent  set  of  assumptions. 
Such  a  code  should  be  able  to  treat  all  currently  known  blade/root  mechanisms  and,  at  the 
same  time,  model  configurations  that  do  not  yet  exist.  It  should  be  possible  to  construct  a 
new  configuration  with  simple  building  blocks  and  with  no  artificial  limitations  on  the  pro¬ 
cess.  For  maximum  flexibility  in  treating  these  different  configurations,  the  finite-element 
method  is  the  preferred  approach.  Moreover,  the  existence  of  many  different  blade/hub 
configurations  for  helicopters  requires  a  capability  to  analyze  arbitrary  configurations  of 
structures,  parts  of  which  may  be  rotating.  Thus,  the  code  should  employ  the  multibody 
philosophy. 


1.2.  Approach 

To  overcome  the  aforementioned  limitations  of  the  existing  methods  of  aeroelastic 
stability  analysis,  the  General  Rotorcraft  Aeromechanical  Stability  Program  has  been  de¬ 
veloped.  GRASP  combines  the  finite-element  and  multibody  approaches,  and  incorporates 
multiple  levels  of  substructures  to  provide  a  powerful  tool  for  rotorcraft  analysis.  The  de¬ 
sign  of  GRASP  is  based  on  the  concept  of  a  collection  of  flexible  and  rigid  bodies  connected 
in  an  arbitrary  manner.  Libraries  of  elements,  constraints,  and  solution  algorithms  appro¬ 
priate  for  the  helicopter  aeroelastic  stability  problem  were  designed  and  built  into  the 
program. 

The  clement  library  promotes  the  modeling  of  the  blades  as  beams;  construction  of 
arbitrary  mechanisms  to  treat  blade/root  kinematics  with  beam  elements  and  rigid  bodies; 
treatment  of  the  fuselage  as  either  a  rigid  body,  a  collection  of  beam  elements,  or  a  modal 
representation  obtained  from  some  other  source;  and  treatment  of  both  static  and  dynamic 
induced  inflow  by  means  of  blade-element/momentum  theory.  The  constraint  library  allows 
arbitrary  connections  between  elements,  includes  constraints  that  allow  for  compliance  in 
the  constrained  relative  motion  between  elements,  and  includes  constraints  that  allow 
the  connection  of  rotating  and  nonrotating  substructures.  None  of  the  constraints  in  the 
library  use  any  kineinatical  approximations,  such  as  small  angle  assumptions.  The  solution 
procedures  include  nonlinear  static  equilibrium  and  linearized  stability  about  equilibriu:  ., 
both  presently  limited  to  the  hovering  flight  condition. 

It  should  be  noted  that  these  physical  modeling  assumptions  and  solution  procedures, 
while  adequate  for  aeromechanical  stability  analysis  in  axial  flight  and  ground  contact, 
are  not  adequate  for  a  comprehensive  rotorcraft  dynamic  analysis  as  defined  by  Johnson 
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(ref.  9).  The  analysis  methodology  used  in  GRASP,  although  a  viable  approach  for  ap¬ 
plication  to  nonlinear  dynamics  in  forward  flight,  would  require  considerable  effort  to  be 
implemented  in  GRASP. 

Several  very  desirable,  but  not  required,  features  of  a  general-purpose  code,  have 
been  incorporated  in  GRASP.  1)  The  accuracy  of  the  analysis  may  be  increased  without 
having  to  add  more  elements.  The  aeroelastic  beam  finite  element  developed  specifically 
for  GRASP  uses  a  variable-order  (or  p-version)  approach,  which  is  based  on  high-order, 
orthonormal,  polynomial  displacement  functions  (refs.  17  and  18).  2)  As  much  as  possible, 
the  equations  of  motion  are  formed  by  the  program  internally,  minimizing  the  possibility 
of  human  error  in  the  equations.  3)  The  user  interface  is  capable  of  handling  a  general 
problem  without  having  to  be  supplied  with  the  form  of  the  equations  of  motion  or  even 
the  number  of  degrees  of  freedom.  4)  Both  large  and  small  problems  can  be  modeled 
with  the  same  code.  Thus,  the  number  of  degrees  of  freedom  is  not  fixed  a  priori.  This 
feature  not  only  requires  a  great  deal  of  flexibility  in  assembling  the  system  equations  of 
motion,  but  also  requires  that  data  be  structured  and  managed  in  core  with  a  flexibility 
not  inherent  in  FORTRAN  (ref.  19). 


2.  SOLUTION  APPROACH 


GRASP  is  specifically  designed  to  provide  a  tool  for  determining  the  equilibrium 
deflections  and  aeroelastic  stability  of  arbitrary  rotorcraft  configurations  in  hover  or  vertical 
flight.  A  GRASP  rotorcraft  model  is  considered  to  be  an  aeroelastic  system  consisting  of 
a  structural  system,  portions  of  which  may  be  rotating  relative  to  one  another,  and  a 
moving  air  mass  with  which  the  structure  interacts.  All  parts  of  the  model  may  be  subject 
to  forces  and  externally  applied  constraints.  The  position  of  any  point  on  the  structure  or 
the  air  velocity  at  any  point  in  the  flow  field  relative  to  an  inertial  frame  of  reference  may 
be  determined  by  solving  a  system  of  partial  differential  and  boundary  value  equations. 
These  equations  are  obtained  from  the  laws  of  fluid  and  structural  mechanics,  and  from 
the  constitutive  properties  of  the  materials  in  the  structure  and  the  air. 

In  vertical  flight,  hover,  or  ground  contact  a  rotorcraft  can  assume  a  steady-state 
equilibrium  configuration  when  the  airflow,  gravity,  and  the  rotor  spin  axis  are  aligned; 
and  when  the  angular  velocity  of  the  rotationally  isotropic  rotor  is  constant.  In  this 
restricted  case  where  the  structure  is  not  subject  to  time-varying  forces,  it  is  possible  to 
eliminate  explicit  time  dependence  from  the  equations.  This  steady-state  equilibrium  can 
be  considered  to  be  static  when  contrasted  to  the  more  general  periodic  equilibrium  found 
in  forward  flight  problems.  The  steady-state  equilibrium  configuration  is  characterized  by 
a  time-invariant  deformation  in  the  nonrotating  portions  of  the  rotorcraft,  a  steady  flow  of 
air  through  the  rotor  disk,  and  time  invariant  deformations  of  the  rotor  blades  with  respect 
to  a  rotating  reference  frame.  The  steady-state  solution  then  calculates  the  equilibrium 
values  of  all  of  the  model  generalized  coordinates  and  generalized  forces. 

The  equations  of  motion  for  the  continuous-structure  portions  of  the  structure  are 
discretized  by  means  of  variable-order,  finite-element  shape  functions.  The  equations  for 
the  structure  then  become  a  system  of  nonlinear,  ordinary  differential  equations.  It  is 
possible,  as  indicated  above,  to  eliminate  all  explicit  dependence  on  time  from  the  equations 
for  the  restricted  case  of  axial  flight  or  ground  contact.  A  linearized  system  of  equations 
may  then  be  calculated  by  taking  small  perturbations  about  the  static  equilibrium  state. 
The  stability  problem  is  defined,  therefore,  by  a  second-order  system  of  linear  equations 
with  constant  coefficients. 

For  infinitesimally  small  perturbations  about  a  previously- calculated,  steady-state 
configuration,  the  dynamic  motion  of  the  rotorcraft  can  be  represented  as  a  linear  combi 
nation  of  complex  eigensolutions.  Since  the  aeroelastic  stability  of  the  rotorcraft  can  be 
determined  directly  from  the  eigenvalues,  the  primary  objective  of  GRASP  can  be  satisfied 
by  computing  these  eigensolutions.  The  frequency  and  damping  information  in  the  eigen 
values  and  the  modal  information  in  the  eigenvectors,  which  can  also  be  obtained  from  the 
eigensolutions,  facilitate  the  user’s  understanding  of  the  dynamics  of  the  rotorcraft. 

The  eigensolution  provides  the  complex  eigenvalues  and  eigenvectors  for  all  model 
degrees  of  freedom  associated  with  the  equations  of  motion  Mq  f  Cq  f  Kq  0  which  have 
been  linearized  about  a  steady  state  deformation.  These  equations  are  often  referred  to 
ns  being  “asymmetric”  because  of  the  nunsymmetry  due  to  aerodynamics  contributions 
to  the  coefficient  matrices  C  and  K.  The  coefficient  matrix  Af,  which  is  both  symmetric 
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and  positive-definite,  contains  contributions  from  the  mass  of  the  structural  model  and 
from  the  “apparent  mass”  of  the  air.  The  coefficient  matrix  C  contains  contributions  from 
structural  and  aerodynamic  damping  and  inertial  forces.  The  coefficient  matrix  K  contains 
contributions  from  structural  stiffness  and  effective  stiffness  from  aerodynamic  and  inertial 
forces.  Like  the  steady-state  solution,  this  solution  requires  that  the  model  correspond  to 
a  physical  system  which  is  not  subject  to  time-varying  forces.  Currently,  the  asymmetric 
cigensolution  must  be  computed  by  using  the  steady-state  solution  obtained  for  an  identical 
model.  This  solution  procedure  prohibits  one,  for  example,  from  obtaining  the  steady-state 
deformations  of  an  isolated  blade,  then  applying  that  solution  to  a  coupled,  rotor/fuselage 
configuration. 
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3.  MODELING  APPROACH 


In  order  to  form  a  mathematical  representation  of  a  structure  that  may  contain  bodies 
which  are  experiencing  large  kinematic  motions  relative  to  one  another,  it  is  necessary  to  be 
able  to  write  the  full,  nonlinear  equations  of  motion  for  the  structure.  The  fundamentals  of 
the  approach  used  in  GRASP  to  derive  these  equations  are  adapted  from  methods  that  were 
originally  developed  for  spacecraft  applications  (ref.  20).  For  the  types  of  structures  that 
GRASP  is  designed  to  represent,  additional  emphasis  ha?  'seen  placed  on  using  multiple 
levels  of  substructures  to  model  a  structure. 

The  first  step  in  modeling  a  structure  in  this  manner  is  to  decompose  the  structure 
(called  the  parent)  into  a  set  of  subordinate  substructures  (called  children),  each  of  which  in 
turn  may  also  be  decomposed  into  a  set  of  child  substructures.  This  decomposition  process 
continues  until  every  substructure  has  been  decomposed  Into  simple  structural  elements. 
The  lowest-level  substructures  (t.e.  those  with  no  children)  are  called  elements.  The  result 
of  this  method  of  modeling  a  structure  is  a  hierarchically-ordered  set  (tree)  of  substructures 
(fig.  1)  that  has  the  complete  structure  at  the  root  and  elements  at  each  of  the  leaves. 
Under  this  modeling  scheme,  a  parent  substructure  may  have  any  number  (including  zero) 
of  child  substructures  but  only  one  parent  substructure.  The  only  substructure  without  a 
parent  is  the  complete  structure,  which  is  at  the  root  of  the  tree. 


Model-type  subsystem:  1 
System-type  subsystems:  2,  3 
Element-type  subsystems:  4,  5, 6,  7 


Figure  1.  Hierarchical  substructure  tree. 


The  hierarchical  model  representation  implemented  in  GRASP  allows  great  generality 
in  the  types  of  configurations  that  can  be  analyzed,  and  permits  essentially  arbitrary  kine 
matic  motions  of  components  relative  to  one  another.  This  general  framework,  along  with 
a  software  design  that  emphasizes  the  use  of  libraries  for  constraints,  elements,  solutions, 
and  so  on,  means  that  the  capabilities  and  limitations  of  the  program  are  those  associated 
with  the  members  of  the  libraries,  not  with  the  program  in  general. 
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3.1.  Subsystems 


In  GRASP,  substructures  are  abstracted  into  subsystems.  Each  substructure  is  then 
represented  by  a  subsystem,  which  may  be  classified  according  to  its  position  in  the  hi¬ 
erarchy  (fig.  2).  The  subsystem  representing  the  complete  structure  (or  model)  is  called 
a  model-type  subsystem.  Substructures  having  no  children  (elements)  are  represented  by 
element-type  subsystems.  The  remaining  subsystems,  those  having  a  parent  and  at  least 
one  child,  are  represented  by  system-type  subsystems.  To  represent  the  substructures  that 
make  up  the  model,  subsystems  serve  several  functions.  First,  they  contain  the  complete 
definitions  of  the  substructures  that  they  represent.  Second,  they  are  repositories  for  the 
generalized  coordinates,  generalized  forces,  and  dynamic  matrices  associated  with  the  sub¬ 
structures.  Finally,  t’  ,y  serve  as  the  basic  units  of  the  hierarchical  organization,  which 
is  an  integral  part  of  the  computational  process  of  transforming  the  parent  generalized 
coordinates  to  the  child  generalized  coordinates,  and  transforming  the  child  generalized 
forces  to  the  parent  generalized  forces. 


Figure  2.  Hierarchical  subsystem  tree. 


Subsystems,  in  general,  may  contain  the  following:  a  frame  of  reference,  a  set  of  nodes, 
a  set  of  generalized  coordinates,  and  a  set  of  constraints.  Each  of  these  entities  perforins 
a  different  function  within  the  subsystem,  and  will  be  described  in  the  following  sections. 
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3.1.1.  Frames  of  Reference 


Every  subsystem  in  a  GRASP  model  (with  the  sole  exception  of  the  air  mass  element) 
has  a  frame  of  reference  associated  with  it.  The  frame  of  reference  is  not  associated  with 
any  material  point  on  the  substructure,  but  instead  serves  as  the  “point  of  view”  for  the 
subsystem.  As  such,  it  establishes  the  coordinate  system  for  that  subsystem.  The  initial 
position  and  orientation  of  a  reference  frame  may  be  selected  to  define  a  coordinate  system 
that  is  natural  for  the  subsystem  ( e.g .,  a  hub-centered  frame  of  reference  might  be  selected 
for  a  subsystem  that  contains  a  helicopter  rotor). 

Since  reference  frames  are  not  physically  connected  to  any  structure,  but  rather  are 
allowed  to  move  freely,  six  degrees  of  freedom  are  associated  with  each  frame.  These 
degrees  of  freedom  define  the  position  and  orientation  of  the  frame  of  reference  for  the 
current  subsystem  relative  to  the  reference  frame  for  its  parent  subsystem. 

In  addition  to  serving  as  a  reference  for  the  subsystem,  the  frame  of  reference  may  be 
used  to  model  the  discrete  motions  of  the  substructure.  This  can  often  lead  to  significant 
simplifications  in  the  equations  of  motion  for  subordinate  subsystems.  For  example,  if  a 
reference  frame  is  attached  to  the  root  of  a  rotating  beam  and  used  to  model  the  rotational 
motion  of  the  beam,  the  equations  of  motion  of  the  beam  itself  need  not  explicitly  include 
the  rotational  motion. 

Since  Newton’s  laws  hold  only  in  an  inertial  reference  frame,  the  model-type  subsystem 
at  the  root  of  the  tree  is  defined  to  be  fixed  in  an  inertial  frame  of  reference.  Therefore, 
while  a  model-type  subsystem  does  have  a  frame  of  reference,  that  reference  frame  has  no 
degrees  of  freedom  associated  with  it  since  it  must  be  inertial.  As  a  result,  the  motions  of 
every  part  of  the  system  can  be  related  to  an  inertial  frame  of  reference. 

3.1.2.  Nodes 

Nodes  arc  used  by  GRASP  to  introduce  degrees  of  freedom  into  a  model.  In  general, 
the  degrees  of  freedom  introduced  by  a  node  may  be  any  generalized  coordinates  that  can 
be  associated  with  a  physically  identifiable  property  of  the  structure.  For  example,  the  set 
of  degrees  of  freedom  for  a  node  could  be  defined  to  be  the  three  rigid-body  translations 
and  the  three  rigid-body  rotations  of  a  point  on  a  structure.  Alternatively,  there  could  be 
a  node  whose  degrees  of  freedom  are  defined  to  be  modal  coordinates. 

Currently,  two  different  types  of  nodes  are  used  by  GRASP:  structural  nodes  and  air 
nodes.  The  structural  nodes  provide  the  measures  for  the  local  displacement  and  rotation 
of  a  structure.  They  move  with  the  deformation  of  the  structure  and  may  be  conceptualized 
as  massless,  infinitesimal,  rigid  bodies  that  are  physically  attached  to  the  structure.  The 
air  nodes  define  the  induced  inflow  velocity  field  through  a  helicopter  rotor.  The  degrees  of 
freedom  for  the  air  node  are  measures  of  the  velocity  distributions  around  the  rotor  disk. 
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3.1.3.  Constraints 


The  constraints  act  as  a  sort  of  “glue”  that  holds  a  model  together.  Constraints 
arc  used  to  model  both  physical  constraints  ( e.g .,  pins,  gimbals,  and  clamps),  and  to 
eliminate  the  dependent  degrees  of  freedom  that  have  been  introduced  into  the  model.  An 
example  of  a  physical  constraint  would  be  the  clamped  boundary  condition  at  one  end  of 
a  cantilever  beam.  That  end  of  the  beam  is  modeled  by  constraining  the  node  at  one  end 
of  the  beam  to  have  no  translational  or  rotational  motion.  Now  consider  two  frames  of 
reference  that  are  defined  to  move  as  if  they  are  rigidly  connected  to  one  another.  For  this 
system,  there  are  twelve  degrees  of  freedom  (six  for  each  frame),  but  only  six  of  them  are 
independent.  Therefore,  a  constraint  must  be  defined  to  remove  the  dependei  t  degrees  of 
freedom.  In  general,  the  set  of  constraints  for  a  subsystem  must  be  sufficient  to  reduce  the 
total  number  of  degrees  of  freedom  to  only  the  independent  degrees  of  freedom  for  that 
subsystem.  Similarly,  for  the  complete  model,  all  dependent  degrees  of  freedom  must  be 
eliminated. 

All  of  the  constraints  implemented  in  GRASP  are  based  on  purely  kinematical  rela¬ 
tionships.  There  are  no  restrictions  to  small  or  moderate  displacements  or  rotations  in  any 
of  the  constraint  equations.  However,  it  is  necessary  to  avoid  the  singularity  that  occurs 
for  deformation-induced  rotations  of  180°  .  This  singularity  arises  as  a  result  of  using 
finite-rotational  kinematics  that  are  based  on  Rodrigues  parameters  (ref.  21). 

The  constraints  in  GRASP  are  implemented  at  two  levels:  the  program  level  and 
the  user  level.  The  constraint  “primitives”  are  found  at  the  program  level.  These  simple 
constraints  provide  a  basic  set  of  connections  among  generalized  coordinates,  frames,  and 
nodes.  At  the  user  level,  these  primitive  constraints  are  combined  to  provide  the  user  with 
physically  meaningful  constraints  between  structural  elements.  For  example,  the  rigid- 
body  mass  connectivity  constraint,  which  is  used  to  attach  a  rigid-body  mass  element  to 
a  structure,  is  a  combination  of  a  primitive  constraint  between  frames  and  a  primitive 
constraint  between  nodes. 

In  order  to  provide  a  full  set  of  constraints,  the  constraint  library  in  GRASP  includes 
several  different  classes  of  constraints.  These  include  constraints  between  two  frames,  con 
straints  between  two  nodes,  constraints  between  generalized  coordinates,  and  constraints 
between  a  frame  and  a  node. 


3.2.  Elements 

Elements  are  subsystems  that  have  no  children.  In  addition  to  frame  and  nodal 
degrees  of  freedom,  they  may  also  have  additional,  non  nodal  generalized  coordinates. 
Computationally,  the  elements  are  the  primary  source  of  virtual  work  in  the  structure. 
For  steady-state  problems,  the  elements  return  the  generalized  forces  associated  with  a 
given  set  of  generalized  displacements.  Fur  perturbation  problems,  the  elements  return 
the  element  coefficient  matrices.  These  matrices  are  determined  from  the  perturbations 
in  generalized  forces  resulting  from  perturbations  in  the  generalized  coordinates  and  their 
time  derivatives. 
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3.2.1.  Rigid-body  Mass 

The  rigid-body  mass  element  represents  a  rigid  body  that  is  subject  only  to  inertial 
and  gravitational  forces.  It  has  a  single  structural  node  that  is  located  at  the  mass  center, 
and  its  axes  are  aligned  with  the  principal  axes  of  the  body.  The  frame  of  reference  for  the 
rigid-body  mass  element  coincides  with  the  nodal  coordinates  in  their  undeformed  state. 

3.2.2.  Air  Mass 

The  air  mass  element  models  the  momentum  air  flow  through  an  axisymmetric  rotor 
disk.  The  degrees  of  freedom  associated  with  this  element  arc  introduced  through  a  single 
air  node.  Since  the  air  mass  element  is  defined  to  be  fixed  in  inertial  space,  the  frame 
degrees  of  freedom  are  suppressed.  For  steady-state  problems,  the  residuals  corresponding 
to  the  uniform  inflow  velocity  and  the  radial  velocity  gradient  are  calculated  from  momen¬ 
tum  considerations  (ref.  22).  For  the  asymmetric  eigenproblem,  only  the  momentum  terms 
(ref.  23)  involving  uniform  and  first-harmonic,  cyclic  perturbations  of  the  inflow  velocity 
contribute  to  the  element  coefficient  matrices. 

3.2.3.  Aeroelastic  Beam 

The  aeroelastic  beam  element  represents  a  slender,  nonuniform  beam  (without  shear 
deformation)  that  is  subject  to  elastic,  inertial,  gravitational,  and  aerodynamic  forces. 
The  primary  assumption  in  the  derivation  of  the  element  equations  (ref.  24)  is  that  strains 
remain  small  relative  to  unity.  There  are  no  small-angle  approximations  made  and  all 
kinematically  nonlinear  effects  are  included.  One  current  limitation  is  that  orientation 
angles  (ref.  21)  (of  type  body-three:  1  2-3)  arc  used  in  the  description  of  finite  rotation 
inside  the  beam  element.  Thus,  rotations  due  to  the  deformation  of  beam  elements  may 
not  exceed  90°  . 

The  aeroelastic  beam  element  degrees  of  freedom  come  from  a  frame  of  reference  that 
coincides  with  the  root  of  the  element  in  its  undeformed  state,  structural  nodes  at  the 
root  and  tip,  an  air  node,  and  a  set  of  internal  degrees  of  freedom.  The  internal  degrees 
of  freedom  result  from  the  higher-order  polynomials  that  may  be  used  to  increase  the 
accuracy  of  the  beam  deformation  calculations.  When  no  internal  degrees  of  freedom  are 
specified,  the  aeroelastic  beam  is  an  Euler- Bernoulli  beam  in  which  the  axial  and  torsional 
deformations  in  excess  of  a  built-in  pretwist  are  represented  by  linear  polynomials,  while 
the  bending  deflcct'ons  are  represented  by  cubic  polynomials.  The  method  of  adding 
internal  degrees  of  freedom  to  improve  the  accuracy  of  an  element  is  more  convenient  than 
adding  elements,  and  is  also  more  efficient  (ref.  17)  given  the  same  number  of  degrees  of 
freedom.  Internal  degrees  of  freedom  may  be  added  selectively  to  reflect  the  dynamics  of 
the  clement.  For  example,  if  a  beam  is  very  stiff  in  bending  and  extension  but  soft  in 
torsion,  additional  torsional  degrees  of  freedom  may  be  added  without  having  to  include 
any  more  bending  or  extensional  degrees  of  freedom. 

The  aerodynamic  forces  on  the  beam  element  are  calculated  from  quasi  steady  strip 
theory  using  lift,  drag,  and  moment  coefficients  that  are  piecewise  continuous  functions  of 
the  angle  of  attack.  Spanwise  scale  factors  for  the  lift,  drag,  and  moment  may  be  specified 
to  allow  for  tip  loss  and  othoi  similar  effects.  The  chord  width,  the  pitch  angle  of  the 


12 


zero-lift, -line,  and  the  offset  of  the  aerodynamic  center  from  the  elastic  axis  may  also  vary 
over  the  length  of  the  element.  The  aeroelastic  beam  element  also  calculates  the  blade- 
element  contributions  to  the  induced  velocity,  which  are  combined  with  the  momentum 
contributions  from  the  air  mass  element  elsewhere. 
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4,  SOLUTION  METHODS 


The  solutions  currently  implemented  in  GRASP  allow  the  user  to  calculate  the  steady- 
state  deformations  of  a  structure  under  load,  and  then  to  solve  for  the  eigenvalues  and 
eigenvectors  of  the  deformed  structure.  In  order  to  obtain  a  valid  eigensolution,  the  steady- 
state  deformations  that  are  used  must  be  such  that  the  structure  is  in  equilibrium. 

4.1.  Steady-State  Solution 

The  equations  for  the  steady-state  equilibrium  of  the  model  are  a  set  of  nonlinear, 
algebraic  equations  of  the  form 

Qi=  i  =  l,...,jv  (4.1-1) 

where  the  Qi  are  the  generalized  forces  (residuals),  the  qi  are  the  generalized  coordinates, 
and  N  is  the  number  of  system  degrees  of  freedom.  These  equations  arc  generated  in¬ 
ternally  by  GRASP  at  the  element  level,  and  automatically  assembled  by  the  constraints, 
which  combine  the  contributions  from  the  finite  elements  into  the  final  set  of  equations.  The 
solution  to  this  set  of  equations  is  obtained  through  the  use  of  the  Levenberg-Marquardt 
algorithm.  This  algorithm  minimizes  the  sum  of  the  squares  of  the  residuals  from  the 
steady-state  equations.  The  implementation  in  GRASP  uses  the  IMSL  (ref.  25)  subrou¬ 
tine  ZXSSQ. 

For  problems  involving  the  aeroelastic  beam  element  with  internal  degrees  of  freedom, 
the  solution  algorithm  is  used  at  two  levels.  First,  it  is  used  in  an  outer  iteration  loop  to 
arrive  at  a  solution  to  the  steady-state  equations  for  the  complete  model  (which  excludes 
the  aeroelastic  beam  internal  degrees  of  freedom).  In  addition,  it  is  used  in  a  separate, 
inner  iteration  loop  to  calculate  the  internal  degrees  of  freedom  for  each  aeroelastic  beam 
element.  A  full  inner  solution  for  each  aeroelastic  beam  is  calculated  for  each  iteration  of 
the  outer  solution. 

In  order  to  arrive  at  a  steady-state  solution,  the  residual  forces  on  the  system  must 
be  calculated,  given  a  deformation  state.  The  algorithm  that  is  used  to  calculate  the 
residuals  for  the  top-level  subsystem  in  the  hierarchical  organization  of  the  model  is  based 
on  a  full-order  tree  traversal  (fig.  3).  When  traversing  down  the  tree  (away  from  the  root 
subsystem),  the  st~..e  vector  for  each  child  subsystem  is  calculated  from  that  of  its  parent. 
Also,  the  inertial  motion  of  the  child  subsystem  reference  frame  is  calculated  from  that  of 
the  parent.  Upon  reaching  an  element,  the  state  vector  for  that  element  and  the  inertial 
motion  for  the  element  frame  are  used  to  calculate  the  element  residuals.  Traversing 
back  up  the  tree  (towards  the  root  subsystem)  the  residuals  from  each  child  subsystem 
arc  transformed  into  its  parent  subsystem  and  added  to  the  parent  residuals.  When  the 
traversal  is  complete,  the  residuals  corresponding  to  each  generalized  coordinate  in  the  root 
subsystem  are  known.  The  complementary  processes  of  calculating  the  state  vectors  and 
assembling  the  residual  vectors  are  accomplished  by  using  the  constraints,  which  define 
the  relationships  among  the  degrees  of  freedom  in  the  parent  and  child  subsystems. 
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3. 

/ 


Calculate  generalized  coordinates  (1, 2, 3, 5, 7, 10) 

Assemble  generalized  forces  (4, 6,  8, 9, 11, 12) 

Forces  calculated  in  subsystems  4,  5,  6,  7 

Figure  3.  Steady-state  solution  full-order  traversal. 

The  solution  methods  available  in  the  current  version  of  GRASP  are  restricted  in  that 
the  same  model  must  be  used  for  both  the  steady-state  and  asymmetric  eigenproblem  solu¬ 
tions.  This  creates  a  problem  for  the  steady-state  solution  algorithm  when  a  configuration 
contains  unconstrained  degrees  of  freedom.  This  can  occur  when  a  model  having  both  ro¬ 
tating  and  nonrotating  components  is  being  analyzed.  For  such  a  configuration,  the  cyclic 
degrees  of  freedom  generated  by  the  rotating  constraints  are  unconstrained.  It  can  also 
occur  in  airborne  configurations,  which  suffer  from  the  same  problem  because  their  body 
degrees  of  freedom  are  unconstrained.  To  alleviate  this  problem,  GRASP  currently  marks 
these  unconstrained  degrees  of  freedom  during  the  building  of  the  model,  and  eliminates 
them  from  the  state  vector  used  in  the  minimization  algorithm. 

4.2.  Asymmetric  Eigenproblem  Solution 

The  system  equations  for  the  asymmetric  eigenproblem  can  be  expressed  in  the  famil¬ 
iar  form 

Mq  -}-  Cq  +  Kq  —  0  (4.2—1 ) 

where  the  q's  are  infinitesimal  perturbations  of  the  generalized  coordinates.  The  algorithm 
used  to  assemble  the  coefficient  matrices  for  the  root  subsystem  is  very  similar  to  that 
used  to  calculate  the  steady-state  residualb  in  that  it  also  is  based  on  a  full  order  tree 
traversal  (fig.  4).  However,  while  traversing  down  the  tree,  no  state  vector  calculations 
are  required.  Upon  reaching  an  element,  the  coefficient  matrices  for  that  element  are 
calculated.  During  the  traversal  back  up  the  tree,  the  constraints  are  used  to  assemble  the 
child  subsystem  matrices  into  the  parent  matrices.  At  the  conclusion  of  the  traversal,  the 
coefficient  matrices  for  the  model  subsystem  are  complete. 
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9  10N 


Assemble  subsystem  matrices  (4, 6,  8,  9, 11, 12) 

Calculate  element  matrices  in  subsystems  4, 5, 8, 7 

Figure  4.  Eigensolution  full-order  traversal. 

The  solution  of  this  set  of  equations  is  begun  by  factoring  matrix  M  using  the  Cholesky 
decomosition  algorithm.  The  GRASP  implementation  uses  subroutine  LUCECP  from  the 
IMSL  (ref.  25)  library.  M  then  becomes 


Introducing  the  transformation 


M  =  LLt 


z  =  LT q 


(4.2-2) 


(4.2-3) 


the  mass  matrix  M  can  be  reduced  to  an  identity  matrix  and  the  system  equations  can  be 
written  as 

A 1  +  L~1CL~rz  +  L~:  KL~tz  =  0  (4.2-4) 


Writing  this  system  of  equations  in  first-order  form 


A  0  t  _  0  A 

0  A  V~  -L-'KL~t  -L~1CL~t  y 


where 


*  =  {!} 


(4.2-5) 


(4.2-6) 


Time  may  be  eliminated  by  the  introduction  of 


(4.2-7) 


10 


which  allows  the  extraction  of  eigenvalues  and  eigenvectors  directly  from  the  matrix  on  the 
right-hand  side  of  equation  (4.2-5).  The  dynamic  matrix  is  balanced,  converted  to  Hes- 
senberg  form,  and  then  the  QR  algorithm  is  used  to  obtain  the  eigensolution.  Finally,  the 
eigenvectors  are  transformed  back  to  the  original  coordinate  system  via  the  transformation 

q *  =  L~tz *  (4.2-8) 

GRASP  uses  subroutine  RG  from  the  NASA/Ames  Cray  library  to  calculate  the 
eigensolutions. 
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5.  COORDINATE  SYSTEMS 


In  GRASP,  many  different  coordinate  systems  are  used  to  mathematically  describe 
the  physical  structure.  To  differentiate  among  them,  each  coordinate  system  is  identified  in 
its  undeformed  state  by  a  capital  ietter  (e.g.,  A).  Depending  on  the  context,  an  identifier 
may  refer  either  to  the  coordinate  system  itself  or  to  a  point  located  at  the  origin  of  the 
coordinate  system.  The  addition  of  a  prime  or  a  double-prime  to  the  identifier  indicates 
that  the  designated  coordinate  system  either  is  in  a  state  of  static  equilibrium  (e.g.,  A1)  or 
is  in  a  dynamically  perturbed  state  (e.g.,  A").  With  these  multiple  coordinate  systems,  it 
is  often  desirable  to  use  several  types  of  mathematical  notation  when  deriving  and  writing 
equations.  Not  only  can  the  form  of  the  equations  be  simplified,  but  also  they  can  be 
made  more  readable.  This  section  is  intended  as  an  introduction  to  the  notation  used  in 
the  sections  where  the  equations  are  actually  derived. 

5.1.  Vectors 

Vectors  play  an  important  role  in  coordinate  system  mathematics.  Associated  with  the 
orthogonal  axes  emanating  from  the  origin  of  every  coordinate  system  is  a  set  of  dextral 
unit  vectors.  These  unit  vectors  are  called  the  base  or  basis  vectors  of  the  coordinate 
system.  In  addition,  vectors  are  used  to  define  variables  such  as  position,  velocity,  and 
acceleration.  Three  types  of  notation  are  used  in  writing  vector  expressions  and  operations: 
vector-dyadic  notation,  index  notation,  and  vector  notation. 

5.1.1.  Vector- Dyadic  Notation 

All  vectors  and  dyadics  used  in  GRASP  are  underlined  (e.g.,  F),  and  all  unit  vectors 
are  identified  by  a  circumflex.  The  difference  between  a  vector  and  a  dyadic  should  always 
be  clear  from  the  context  of  its  usage.  For  a  coordinate  system  A,  the  basis  vectors  are 

written  as  6,  ,  where  i  =  1,2,3.  Any  unit  vector  other  than  a  basis  vector  is  denoted  by  e, 
and  may  appear  either  with  or  without  superscripts. 

When  kinematical  quantities  have  coordinate  systems  associated  with  them,  the  rela¬ 
tionship  is  defined  by  using  the  appropriate  superscripts.  For  example, 

Rda  =  position  of  the  origin  of  coordinate  system  B 

with  respect  to  the  origin  of  coordinate  system  A 

V_BA  A  velocity  of  the  origin  of  coordinate  system  B 
with  respect  to  coordinate  system  A 

Aba  =  acceleration  of  the  origin  of  coordinate  system  B 
with  respect  to  coordinate  system  A 

=  angular  velocity  of  coordinate  system  B 
with  respect  to  coordinate  system  A 
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Forces  and  moments  are  significant  in  their  point  of  application  as  well  as  their  source. 
The  notation  adopted  herein  is 

Fa  =force  at  A 
Ma  =moinent  at  A 

For  example,  a  force  and  moment  at  A  contribute  to  a  moment  at  D  according  to  the 
relationship 

Mb=Ma  +  RadxAa  (5. 1.1-1) 

5.1.2.  Index  Notation  for  Vectors 
A  vector  V_  in  the  .4  basis  may  always  be  expressed  as 

V  =  VAib?  (5. 1.2-1) 

where  the  summation  convention  adopted  is  that  repeated  indices  are  always  summed 
over  their  range.  Unless  otherwise  specified,  Latin  indices  assume  the  values  1,2,3;  Greek 
indices  assume  the  values  1,2.  The  subscript  A  in  VAl  indicates  that  the  measure  numbers 
VAi  are  defined  in  the  A  basis. 

Two  symbols  frequently  encountered  in  vector  operations  that  use  index  notation  are 
the  Kronecker  delta  Sij  and  the  Levi-Civita  epsilon  e,yfc  where 


i  £  j 
i  =  j 


(5.1. 2-2) 


{0  any  index  repeated 

+1  cyclic  permutation 

-1  acyclic  permutation 


(5. 1.2-3) 


The  Kronecker  delta  consists  of  the  components  of  the  identity  tensor  in  a  Cartesian 
coordinate  system,  while  the  Levi  Civita  epsilon  consists  of  components  of  the  permutation 
tensor  in  a  Cartesian  coordinate  system.  Some  useful  identities  regarding  both  of  these 
symbols  may  be  found  in  reference  26. 


5.1.3.  Matrix  Notation  for  Vectors 


Using  index  notation,  a  vector  V  may  be  expressed  in  the  A  basis  as  shown  in  equa¬ 
tion  (5.1.2  1).  Since  the  basis  is  identified  by  the  subscript  A,  the  measure  numbers 
themselves  may  be  viewed  as  a  complete  description  of  the  vector.  Thus,  the  column 
matrix  U,i  can  be  defined  to  be 

(  ^  1 

Va  -  <  VM  (5. 1.3-1) 
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as  an  alternate  way  of  expressing  the  vector  V.  The  dot  product  U_-  V_  may  then  be  written 
as 

(  Kii) 

UjVA  =  [UM  UA2  U A 3J  {  V A2 

1^3. 

The  cross  product  of  two  vectors  17  and  V_  may  be  written  as 


(5. 1.3-2) 


u  X  V  =  UAit t  X  vAib? 


- A 

—  £\jk  UAjVAkbi 


(5. 1.3-3) 


=  UAijVAjb? 


This  equation  implies  that  the  measure  numbers  of  the  cross  product  in  the  A  basis  are 
simply  the  elements  of  the  matrix  product  UaVa  where 


(  )„  = 

-tij  k{ 

)fc 

For  example, 

0 

-Ua2 

Ua2  ' 

Ua  = 

Ua  3 

0 

~UAi 

-UA2 

Uai 

0 

(5. 1.3-4) 


(5. 1.3-5) 


There  are  also  several  useful  identities  that  can  be  derived  for  two  column  matrices  a 
and  b 

-T 

a  =  —  a 

ab  =  —  ba 


a1b  =  —  bTa  =  (ab)1 
ab  =  -  a1  bA  +  ba1 
ab  =baT  —  abT 


(5. 1.3-6) 


ab  —  ab  =ba 
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5.2.  Finite  Rotations 


In  many  kinematic  analyses,  rotations  are  assumed  to  be  either  infinitesimal  or  moder¬ 
ate  in  size.  These  assumptions  allow  certain  simplifications  in  the  kinematical  relationships, 
but  constrain  the  range  of  applicability  of  the  analysis.  In  GRASP,  no  such  assumptions 
are  made  and  all  rotations  are  assumed  to  be  of  arbitrary  size  (finite).  Finite  rotations  are 
expressed  in  four  ways  in  GRASP: 

(1)  direction  cosines, 

(2)  Euler  rotations, 

(3)  Tait-Bryan  orientation  angles,  and 

(4)  Erder-Rodrigues  parameters.  Internally,  GRASP  expresses  all  finite  rotations  in  terms 
of  direction  cosine  matrices.  For  the  convenience  of  the  user,  any  of  the  other  three 
methods  may  be  used  to  specify  the  input  to  GRASP.  Since  there  are  significant 
differences  in  the  algorithms  used  to  compute  the  direction  cosine  matrix,  all  three  of 
the  other  representations  are  also  discussed  in  detail. 


5.2.1.  Direction  Cosines 

When  a  coordinate  system  B  undergoes  an  arbitrary  rotation  relative  to  coordinate 
system  A,  the  basis  vectors  are  related  by  the  equation 


(5. 2.1-1) 


where  the  superscripts  arc  coordinate  system  identifiers,  not  indices.  The  matrix  of  direc¬ 
tion  cosines  CBA  is  orthonormal  such  that 


CBACAB  -  CABCBA  =  A  (5.2. 1-2) 

It  should  be  noted  that  the  form  of  the  matrix  of  direction  cosines  used  in  this  manual  is 
the  transpose  of  that  developed  in  reference  21. 

Similarly,  with  this  notation  it  is  easy  to  show  that  a  basis  change  for  any  kinematical 
vector  can  be  performed  by  changing  the  subscript  and  multiplying  by  the  matrix  of 
direction  cosines  for  the  bases. 

VB  =  CbaVa  (5.2.1-3) 

Note  that  for  kinematical  vectors  the  superscripts  are  unaffected  by  these  operations. 


5.2.2.  Euler  Rotations 

If  coordinate  system  B ,  initially  coincident  with  A,  rotates  about  a  unit  vector  e  fixed 
in  A  by  an  angle  0  (fig.  5)  then  the  matrix  of  direction  cosines  can  be  written  as 


CBA  =  A  cos#  +  Cjxca1  (1  -  cos  0)  —  eA  sin0 


Tt  i 


(5.2.2-f ) 


where 


(5. 2. 2-2) 
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Figure  5.  Euler  rotation. 

5.2.3.  Tait-Bryan  Orientation  Angles 

Consider  two  coordinate  systems  A  and  B  with  coincident  basis  vectors  bf  and  bf . 
Let  the  orientation  of  B  with  respect  to  A  change  as  follows  (fig.  6): 

(1)  Perform  an  Euler  rotation  of  B  about  e  =  ( j  —  1,2,  or  3)  by  an  angle  Qj\ 

. D 

(2)  Perform  an  Euler  rotation  of  B  about  e  =  bk  (k  =  1, 2,  or  3,  k  ^  j)  by  an  angle  0*; 

« B 

(3)  Perform  an  Euler  rotation  of  B  about  e  =  bt  (l  =  1,2,  or  3,  /  /  kyl  ^  j)  by  an  angle 
0f. 

The  final  orientation  of  B  relative  to  A  depends  both  on  the  magnitudes  of  0j,  92 ,  and  93 
and  the  sequence  j-k-l.  Details  of  this  type  of  transformation  may  be  found  in  reference  21 
where  Tait-Bryan  angles  are  classified  as  orientation  angles  of  type  body-three.  For  the 
rotation  sequence  1-2-3  the  matrix  of  direction  cosines  is  calculated  as  follows: 

C3  S3  0  C2  0  —32  10  0 

CDA  -  -s3  c3  0  0  1  0  0  cj  sj 

0  0  1  s2  0  c2  0  —si  cj 

(5. 2. 3-1) 

C2C3  C]S3  +  SiS2C3  SiS3  —  CiS2C3 

=  —  C2s3  Cj  C3  —  Sj  S2  S3  S]C3  —  C]S2S3 
S2  — ^1 C3  C1C2 


v  tere 


c,  —  cos0, 
3,  -  sin0, 


(5.2.3  2) 
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Figure  6.  Tait-Bryan  orientation  angles  (1-2-3). 


5.2.4.  Euler- Rodrigues  Parameters 

For  two  coordinate  systems  >1  and  B ,  three  parameters  <£,  =  2e,i,  tan(|)  may  be  used 
to  describe  a  change  in  orientation  (ref.  21).  The  values  of  herein  are  scaled  by  a  factor 
of  2  relative  to  the  Rodrigues  parameters  presented  in  reference  21,  so  that  for  infinitesimal 

values  of  if)l  -  the  rotation  may  be  regarded  as  a  vector  <£,6,  =  <£,&,  with  CBA  -  &  -  <f>. 
The  matrix  of  direction  cosines  is  then  simply 


yDA 


(1  -  ffi)A  -1  -  j 


1  +  ^ 


(5. 2.4-1) 


The  angular  velocity  of  B  relative  /l,  expressed  in  the  B  basis,  can  be  written  as 

(A 


A 

**» 


1  +  ^ 


(5. 2.4-2) 


These  relations  contain  no  trigonometric  functions  and  are  easilj  expressed  in  a  shorthand 
matrix  notation.  Furthermore,  a  simple  inverse  transformation  exists  s  that  given  CBA, 
the  values  of  <j)  may  be  obtained  from 


DA 


$ » 


i  •' 'eg* 


(5.2.4-3) 


where  CBA  is  the  trace  of  CBA.  Given  and  SlBA,  <i>  can  be  obtained  from 


</>  (A 


2 


WS 


BA 


(5.2.4  4) 


Transformations  between  Euler  Rodrigues  parameters  arrd  direction  vusines  (or  angular 
rates)  are  ver\  simple  relative  t* •  the  transformations  required  for  Tait  Bryan  angles. 
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5.3.  Angular  Velocity  and  Virtual  Rotation 

The  measure  numbers  for  the  angular  velocity  of  coordinate  system  B  relative  to 
coordinate  system  A  expressed  in  the  A  basis,  ilBA,  may  be  determined  from  the  addition 
theorem  discussed  in  reference  27.  They  can  be  related  to  the  time  derivative  of  the  matrix 
of  direction  cosines  as  follows: 

CBA  =  ■ ACBA  =  -CBASl%A  (5.3-1) 

By  virtue  of  the  Kirchhoff  kinetic  analogy  (ref.  28),  CBA  in  Eq.  (5.3-1)  may  be  replaced 
with  SCBA,  and  f lBA  with  Sij)BA.  The  expression  for  the  components  of  virtual  rotation 
of  B  in  A  then  becomes 

6Cba  =  -6i>BBACBA  =  -Cba<4BaA  (5.3-2) 


The  corresponding  virtual  rotation  vector  SipBA  is  used  in  determining  the  virtual 
work  due  to  applied  moments.  The  components  of  virtual  rotation  may  be  obtained  from 
any  expression  involving  the  angular  velocity  in  a  manner  identical  to  that  used  to  obtain 
equation  (5.3-2)  from  equation  (5.3-1). 

Similarly,  infinitesimal  perturbations  of  the  rotation  vector  can  be  obtained  by  sub¬ 
stituting  CUA  for  CDA  and  8% A  for  in  equation  (5.3-1). 

CBA  =  -0BACDA  =  ~Cba6ba  (5.3-3) 


5.4.  Velocity,  Acceleration,  and  Virtual  Displacement 

Velocity  and  acceleration  vectors  are  obtained  by  applying  the  superposition  theo¬ 
rems  discussed  in  reference  27.  The  calculation  of  the  velocity  and  acceleration  vectors 
is  fundamentally  nothing  more  than  the  differentiation  with  respect  to  time  of  a  position 
vector  in  (he.,  relative  to)  some  coordinate  system.  It  is  often  necessary  to  determine  the 
time  derivative  of  a  vector  in  coordinate  system  B,  when  the  derivative  is  known  only  in 
coordinate  system  A.  Given  an  arbitrary  position  vector  R  and  its  first  and  second  time 
derivatives  in  A,  the  first  and  second  time  derivatives  of  R  in  B  may  be  determined  from 
the  following  expressions. 

°t£=  at&  '-Qab  x  R 

at  at 


B 


"t  (  dt&) 

dt  \  dt  ) 

A  f  (  Aix  +  nAB  xr) 
dt  V  dt~  ~  -) 

,  d2  .  d 

17i$+  xR  + 

dt 1  dt 


+  qab  x  ^  A-~R  +  nAB  xrJ 

2Qab  x  A  4r  +  nAB  x  (Qab  x  R) 
dt 


(5.4-1) 
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The  Kirchhoff  kinetic  analogy  (ref.  28),  can  also  be  applied  to  equation  (5.4-1)  to 
obtain  the  virtual  displacement  vector.  Time  derivatives  in  B,  Bji(  )>  are  replaced  with 
B5(  );  velocity  vectors  in  A,  Aj[(  ),  are  replaced  with  virtual  displacments  in  A,  AS(  ); 
and  angular  velocity  vectors  £Lab  are  replaced  with  virtual  rotations  S^AB . 

b6R  =  ASR  +  8±ab  x  R  (5.4-2) 
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6.  SUBSYSTEMS 


As  described  in  Section  3,  the  physical  structure  that  is  being  modeled  by  GRASP  is 
broken  down  into  a  hierarchy  of  substructures.  Each  of  those  substructures  is  represented 
in  GRASP  as  a  subsystem.  Every  subsystem  in  the  model  is  in  turn  composed  of  a  set  of 
components  which  may  include  a  frame  of  reference,  a  set  of  nodes,  a  set  of  constraints, 
and  a  set  of  child  subsystems.  It  is  the  interrelationship  among  these  components  that 
allows  the  construction  of  the  equations  of  motion  for  each  subsystem. 

6.1.  Frames  of  Reference 

The  position  of  the  frame  of  reference  F  for  a  child  subsystem  relative  to  the  frame 
of  reference  5  for  the  parent  subsystem  is  defined  as  RFS ,  and  the  orientation  (direction 
cosines)  of  the  child  subsystem’s  frame  relative  to  the  parent’s  frame  is  defined  as  CFS 
(fig.  7).  Since  Newton’s  laws  apply  only  in  inertial  frames  of  reference,  all  equations  of 
motions  must  be  written  relative  to  an  inertial  reference  frame.  Therefore,  it  is  essential 
to  have  a  method  of  transforming  back  to  the  inertial  frame  from  any  subsystem  frame  in 
the  model.  If  the  position  and  orientation  of  the  parent’s  reference  frame  5  are  defined 
relative  to  an  inertial  reference  frame  7,  the  inertial  position  and  orientation  of  any  child’s 
reference  frame.  F  can  be  determined  from  the  parent’s  reference  frame  S  by  applying  the 
following  equations  recursively. 


RFI  =RrS  -|  RSI 
CFI  =CFSCSI 


(6.1-1) 


In  addition  to  the  position  and  orientation  of  any  reference  frame  relative  to  the  inertial 
reference  frame,  it  is  necessary  to  know  the  inertial  motion  of  every  subsystem  frame.  A 
subsystem  reference  frame  may  experience  accelerations  relative  to  the  inertial  frame  if  it 
or  any  of  its  direct  ancestors  is  experiencing  translational  accelerations  or  rotation  motions. 
Thus,  if  the  inertial  motion  of  the  parent’s  frame  of  reference  5  is  known,  then  the  velocity, 
angular  velocity,  and  acceleration  of  the  child’s  frame  F  can  be  obtained  from  the  following 
equations: 

vFI  =zFS  +  vs  1  +  nS!  x  rfs 

Qf!  =Q.si  +  Qfs  (6.1-2) 

AFl  =Afs  +  A5'  X  (QSI  X  Rfs) 
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Figure  7.  Frames  of  reference. 

When  expressed  in  the  appropriate  bases,  these  equations  become  (in  matrix  form) 

vp  =CFS(VFS  +  vg1  +  nf'Rf5) 

o.£-r  =c,F5(nf/  +  of5)  (6.1-: 

afi  =cfs(affs  +  Ai1  +  nf'nfRf 5) 


Note  that  in  the  current  version  of  GRASP,  it  is  assumed  that  YFS  —  AFS  —  0. 


Frames  also  possess  six  rigid-body  degrees  of  freedom.  Thus,  while  frames  are  not 
physically  attached  to  the  structure,  they  may  move  relative  to  one  another  in  space.  In  the 
case  of  steady-state  deformations,  these  six  degrees  of  freedom  include  three  translations 
along  the  deformed  frame  basis  vectors  and  the  three  Euler-Rodrigues  parameters  for 
angular  displacements.  The  steady-state  displacement  vector  for  frame  F  is 


Rf'f  ---  R?fk 


(0.1-4) 


The  steady-state  frame  rotations  are  expressed  in  terms  of  <f)['F ,  and  the  direction 
cosines  of  the  deformed  frame  coordinate  axes  F'  with  respect  to  the  undeformed  coordi 
nate  axes  F  arc  written  as  CFj  F.  In  matrix  notation,  the  steady-state  frame  state  vector 
is 


qp> 


-  [  1 
“  l  if”  / 


(6.1-5) 


For  dynamic  perturbations  about  the  steady-state  condition,  the  displacement  vector 


is 


rf"f' = if" 


(6.1-6) 


The  dynamic  perturbations  of  the  frame  rotations  arc  expressed  in  terms  of  infinitesimal 
rotations  Op„f  ,  for  which  the  direction  cosines  (ref.  21)  are 


zF"F' 

—  A  —  8  pn  j 


(6.1-7) 


In  matrix  notation,  the  dynamic  perturbation  frame  state  vector  is  then 

_ / *£f \ 

qr"  ~  \  if,'.?  I 


(6.1-8) 


The  virtual  displacements  for  the  steady  state  and  dynamic  formulations  are  simply 
variations  of  the  displacement  coordinates 


snF'F  Mf;Fi? 

SRf"f'  =6  RpfC 

and  the  virtual  rotations  are  variations  of  the  rotational  degrees  of  freedom 

,  i.-'  ,  e"  c-' ;  F" 


(6.1-9) 


(6.1  10) 
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8.2.  Nodes 


Nodes  arc  used  by  GRASP  to  model  the  kinematics  of  a  structure,  and  their  degrees 
of  freedom  are  representative  of  the  physical  states  of  that  structure.  The  position  and 
orientation  of  any  node  is  defined  relative  to  the  frame  of  reference  for  the  subsystem  in 
which  the  node  resides.  Thus,  for  a  node  JV  in  a  subsystem  with  reference  frame  F,  the 
position  and  orientation  of  N  with  respect  to  F  arc  RNF  and  CNF  %  respectively.  Two  types 
of  nodes  are  currently  used  in  GRASP:  structural  nodes  and  air  nodes  The  kinematics  of 
these  nudes  are  described  in  the  following  sections. 


6.2.1.  Structural  Nodes 

A  structural  node  represents  a  specified  material  point  on  a  structure.  Since  the 
material  point  may  have  up  to  six  degrees  of  freedom,  the  structural  node  also  has  six 
degrees  of  freedom.  For  the  case  of  steady-state  deformations,  these  six  degrees  of  freedom 

«AT 

include  three  translations  along  the  undeformed  nodal  basis  vectors  bx  and  three  Euler- 
Rodrigues  parameters  f> ^  N  for  angular  displacements.  The  nodal  displacement  vector  for 
node  N  is  then 

Rn‘n  =  R^b-  (6.2.1— 1) 

The  direction  cosines  of  the  deformed  nodal  coordinate  axes  N'  relative  to  the  undeformed 
axes  N  are  expressed  as  C ^  N .  Then,  in  matrix  notation,  the  nodal  state  vector  is 

f  PN'.N  1 

W  =:  |  |  (6.2.1-2) 


Note  that  the  nodal  steady- state  degrees  of  freedom  are  referenced  to  the  undeformed  nodal 
basis,  whereas  the  frame  steady-state  degrees  of  freedom  are  referenced  to  the  deformed 
frame  basis. 


For  dynamic  perturbations  about  the  steady-state  condition,  the  displacement  vector 


is 


rk"n'  =  rk;n'i? 


(6. 2. 1-3) 


The  dynamic  perturbations  of  the  nodal  rotations  are  expressed  in  terms  of  infinitesimal 
rotations  0$£  N  ,  for  which  the  direction  cosines  are 


qN"N' 


-N"N' 
=  A  -  0Ni 


(6.2. 1-4) 


In  matrix  notation,  the  dynamic  perturbation  nodal  state  vector  is  then 


qNn  - 


(6.2. 1-5) 


Note  that  the  nodal  dynamic  degrees  of  freedom  are  referenced  to  the  undeformed  nodal 
basis,  whereas  the  frame  dynamic  degrees  of  freedom  are  referenced  to  the  dynamically 
perturbed  frame  basis. 
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The  virtual  displacements  for  the  steady-state  and  dynamic  formulations  are  simply 
variations  of  the  displacement  coordinates 


SRN‘N  =6RH,':Nb" 


Ni  2* 

ii  wiAf 


(6.2. 1-6) 


8Rn"n'  =6rx;n'£ 
and  the  virtual  rotations  arc  variations  of  the  rotational  degrees  of  freedom 


Sri 


N"N‘ 


(6.2.1-7) 


=Srj) 


Ni 


6.2.2.  Air  Nodes 

The  generalized  coordinates  representing  the  axisyinmetric  flowfield  associated  with 
a  helicopter  rotor  are  introduced  into  GRASP  by  means  of  the  air  node.  The  generalized 
coordinates  arc  defined  relative  to  an  inertial  frame  of  reference  I,  and  determine  the 
inertial  air  velocity  at  a  point  Q  as 


UQ'  =  -(Uf  +  nt  +  RqJ  it,  +  ^sV.lK1  (6.2.2-1) 


where  is  an  inertially  fixed  unit  vector  and  A  is  a  coordinate  system  whose  origin  is 
located  at  the  center  of  the  axisymmetric  flowfield.  The  distance  from  the  center  of  flow  r 
can  be  calculated  from 


(6. 2. 2-2) 


For  the  case  of  steady-state  inflow,  (7,  and  7^  represent  the  uniform  inflow  velocity 
and  the  radial  velocity  gradient,  respectively.  The  other  two  generalized  coordinates  have 
no  physical  meaning  under  these  conditions,  and  therefore  are  not  used.  The  air  node  state 
vector  for  steady-state  inflow  is  then 


QA'  = 


(6.2. 2-3) 


To  model  dynamic  inflow,  generalized  coordinates  U 7(2,  and  7(3  represent  the 
collective  and  two  cyclic  velocity  perturbations.  The  air  node  state  vector  for  dynamic 
inflow  is  then 


9  A"  — 


Uf  ' 

Tils  > 


(6. 2. 2-4) 
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7.  CONSTRAINTS 


The  purpose  of  a  constraint  is  to  create  a  dependency  among  generalized  coordinates. 
In  GRASP,  the  dependencies  among  the  generalized  coordinates  are  used  to  eliminate 
dependent  generalized  coordinates  in  favor  of  independent  generalized  coordinates.  In 
the  following  sections,  the  general  formulation  of  a  primitive  constraint  will  be  presented, 
followed  by  the  specific  applications  in  GRASP.  Then,  the  composite  constraints  that  have 
been  constructed  from  the  primitive  constraints  will  be  discussed. 


7.1.  Primitive  Constraints 


Consider  a  set  of  generalized  coordinates  that  are  related  to  one  another  through  a 
constraint.  The  constraint  relationship  g  may  be  written  in  the  special  form 

9e(  =  9i  ( Qrl  i  •  •  •  i  QrNr  )>  (i  =  1, . . . ,  N  )  (7.1—1) 

Thus,  the  generalized  coordinates  related  by  the  constraint  can  be  partitioned  into  two 
sets:  a  set  to  be  eliminated,  ge,  and  a  set  to  be  retained,  qr.  Using  the  constraint  relation¬ 
ship,  the  set  to  be  eliminated  can  be  obtained  directly  from  the  constraint  functions  which 
depend  only  on  the  set  to  be  retained. 

The  virtual  work  for  the  generalized  coordinates  associated  with  the  constraint  is 

Ne  Nr 

IW  =  Y,  * <?.,  +  £  ««  ,,  Qr,  (7.1-2) 

;=i  ;=i 


The  sum  of  the  generalized  forces  Q  associated  with  a  generalized  coordinate  may  differ 
from  zero  for  two  reasons.  First,  during  the  process  of  seeking  an  equilibrium  solution, 
equilibrium  may  not  always  be  satisfied.  In  this  case,  the  sum  of  the  generalized  forces  is 
residual  force  that  is  a  measure  of  the  error  in  the  approximate  solution.  Second,  even  if 
the  complete  system  is  in  equilibrium,  individual  subsystems  may  not  be  in  equilibrium. 
The  generalized  forces  for  these  subsystems  will  be  nonzero. 

Taking  the  variation  of  equation  (7.1-1) 


Nr 


ht,  =  §jr6qn »  (*  —  l,...,JVe) 

i=X  0<lri 


}= 1 


(7.1-3) 


and  substituting  Eq.  (7.1-3)  into  Eq.  (7.1-2) 


Nr  (  Ne  a 

;=i  \  t=i  ir; 


(7.1-4) 


This  relationship  is  used  by  GRASP  to  incorporate  the  contributions  of  the  generalized 
forces  associated  with  the  eliminated  generalized  coordinates  into  the  retained  generalized 
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forces.  During  calculation  of  steady-state  residuals,  the  residuals  associated  with  the  elim¬ 
inated  generalized  coordinates  are  transformed  and  added  to  the  appropriate  residuals  in 
the  parent  substructure’s  residual  vector. 

The  treatment  of  constraints  for  small  perturbations  about  an  equilibrium  state  is 
a  little  more  involved.  For  this  problem,  each  generalized  coordinate  is  assumed  to  be 
the  sum  of  an  equilibrium  value  and  an  infinitesimal  perturbation  from  that  value  (i.e., 
q  =  q  - 1-  q).  Equations  (7.1-1),  (7.1-3)  and  the  generalized  forces  Q  can  all  be  expanded 
in  Taylor  series  about  the  equilibrium  value.  Noting  that  equation  (7.1-1)  is  valid  when 
q  —  q,  expansion  of  equation  (7.1-1)  yields 


Nr  •£— 

Em  .  , 

— + 


(i  =  1 . N') 


Expansion  of  equation  (7.1-3)  yields 


(7.1-5) 


(7.1-6) 


Expansion  of  the  generalized  force,  Q,  for  both  eliminated  and  retained  terms  yields 

Ne  Nr 

Qei  =Qei  -b  e,- 9e,-  +  ^  Leir.qr. , 

j=l  j= 1 

N *  Nr 

Qr{  =Qn  +  ^^Lr. e;-3e;-  +  Lrirjqrj , 
j=i  j=i 

where  the  linear  operator,  L,  contains  the  terms  normally  associated  with  the  mass, 
damping,  and  stiffness  matrices,  -M  jp  -  C  jt  -  K.  Note  that  the  minus  signs  are  present 
in  the  definition  of  L  because  the  generalized  force  is  generally  regarded  as  positive  on 
the  right-hand  side  of  the  dynamical  equation,  whereas  the  linear  coefficient  matrices  are 
regarded  as  positive  on  the  left-hand  side. 

GRASP  calculates  the  M,  C,  and  K  matrices  for  a  subsystem  by  adding  the  contri 
butions  of  each  of  its  children.  The  rows  and  columns  of  the  child  subsystem’s  matrices 
correspond  to  all  of  the  generalized  coordinates  of  the  child.  The  constraints  are  used  to 
eliminate  dependent  generalized  coordinates,  resulting  in  matrices  whose  rows  and  columns 
correspond  to  only  the  retained  generalized  coordinates  of  the  cliild.  The  matrices  elements 
are  then  added  to  the  elements  of  the  parent's  matrii.cs  that  correspond  to  the  child’s  inde 
pendent  degrees  of  freedom.  The  required  transformations  can  be  found  using  the  virtual 
work  for  the  subsystem. 


(»=!,...,*•) 


(7.1-7) 


(*  =  1 . N") 
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An  expression  for  the  virtual  work  from  small  perturbations  about  the  equilibrium 
state  may  be  obtained  by  substituting  equation  (7.1-6)  and  the  eliminated  and  retained 
subsets  from  equations  (7.1-7)  into  the  virtual  work  expression  in  equation  (7.1-2). 


W-E5>,  (£+£53^*.+-)  (e,,+£w,.+Ei.,,*, 


+  y  J  Hn  I  Qn  +  ^  Lr. r.qr.  | 

i=l  \  J=1  j= 1  / 

(7.1-8) 

After  discarding  terms  of  second  or  higher  order,  the  expression  for  virtual  work 
consists  of  a  constant  part  and  two  first-order  parts  in  q.  The  constant  part  is  the  same 
as  equation  (7.1-4),  except  that  it  is  evaluated  for  the  equilibrium  state. 


N'  N'  ,  -5— 

££^,(^..+q,. 

t— I  J— 1  \°qri 


(7.1-9) 


The  first  linear  portion  of  the  virtual  work  is  the  single  term 


NT  N ' 


NT  Nr 


j-  t  k=  1 


(7.1-10) 


The  matrix  K°  represents  the  geometric  stiffness  associated  with  the  constraint.  During 
assembly  of  the  matrices  for  the  parent  substructure,  GRASP  calculates  this  geometric 
stiffness  and  adds  it  to  the  stiffness  matrix  in  the  parent  substructure.  This  extremely 
important  term  is  often  u\erlooked.  For  instance,  a  pendulum,  modeled  as  a  rigid  body 
mass  constrained  to  rotate  about  an  offset  axis  (using  a  screw  constraint)  derives  all  of  its 
stiffness  from  this  geometric  stiffness  term. 

The  remainder  of  the  linear  terms  are 


A"  /V' 


EE^,  (v.  + 


A"  Nr 


EEkI^^E^ 

j  - 1  Jt—  1  \  «=i  Hr>  / 


(7.1-11) 
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After  substituting  equation  (7.1-5)  into  equation  (7.1-11)  for  the  eliminated  perturbation 
coordinates  these  terms  become 


Nr  Nr  (  N‘ 

j=  1  k=  1  V  t=l  ^ri 


N'  -z—  N'  N*  -z-  -s—  \ 

Ey  vgi  .  STST  °9i  7  d9l  \  . 


(7.1-12) 


The  quantity  within  the  parentheses  in  equation  (7.1  12)  can  be  thought  of  as  defining 
a  new  set  of  M,  C ,  and  K  matrices  in  terms  of  the  retained  and  eliminated  portions  of 
the.  original  matrices.  GRASP  calculates  the  new  matrices  and  adds  their  elements  to  the 
elements  of  the  parent  substructure’s  matrices. 

The  definition  of  a  constraint,  follows  from  the  specification  of  the  function  //.  To  obtain 
a  solution  for  a  system  in  equilibrium,  the  matrix  ^  must  be  known.  A  perturbation 

solution,  however,  requires  both  the  matrix  ag9  and  the  geometric  stiffness  matrix  KG .  In 
the  following  constraint  derivations,  matrix  will  be  denoted  by  'R.. 

7.1.1.  Fixed  Frame 

The  fixed  frame  constraint  describes  a  rigid  connection  between  two  frames  of  refer 
ence,  F  and  S.  Regardless  of  the  changes  in  position  and  orientation  relative  to  inertial 
space  that  they  may  undergo,  their  position  and  orientation  relative  to  one  another  re 
mains  constant.  The  current  (child)  frame  F  will  have  its  degrees  of  freedom  eliminated, 
while  the  degrees  of  freedom  for  the  superurdinatc  (parent  )  frame  5  will  be  retained.  In 
GRASP,  this  constraint  is  available  through  the  user  interface. 

Steady- State.  Consider  two  frames  in  their  undeformed  (S  and  F)  states  and  in  their 
steady-state  (S'  and  F1)  configurations  (fig.  8).  The  degrees  of  freedom  of  F  (and  F')  are 
considered  to  be  dependent,  while  those  of  5  (and  5')  are  independent.  The  frames  are 
assumed  to  be  connected  such  that 


where 


Rsf  +  Rff'  f  rf's' 


rf:s'  -  rfs 


Thus, 


£-»  S'  .S'  S  !•'  F  F'  P  1  S' 


rF'F  .  rF'S'  ,  RS'S  rfs 


Rli‘.1  Cn'(J?^'  ♦  Up’)  Cl  vcs  srCs 
CFF  p F’S'pS‘  S^S F 


(7.1.1  1) 


(7.1.1  -2) 


(7.1.1  3) 


(7  1  1  4) 
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Figure  8.  Fixed  frame  constraint. 

Consider  the  virtual  work  performed  in  the  F'  frame  at  F' . 

6W  =  (SRp,F)T  Fp,  +  («*£F)rjtf£  (7. 1.1-5) 

The  virtual  displacements  and  rotations  at  F  are  related  to  those  at  5  such  that 

XT)F'F  _  nF'S'fcnS'S  KnS‘Sr,FS\ 


6Ry  =  C*  b  {6RSS,S  -  SC *  bRs) 


=  cf's'(6R*s  +  S4>ss,scs'srP) 


-  Cf's\6R%s  -  R%?6il>f.s) 


s*f:f  -  Cf's'6^ss:s 

so  that  the  virtual  work  performed  in  the  F'  frame  at  S'  becomes 


(7. 1.1-6) 


pi  c  i  pi  pi  pi 

6W-(SR§,s)  Cs  F  Ff,  + 


(7.1. 1-7) 
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Equations  (7.1. 1-7)  show  the  contributions  of  the  force  and  moment  acting  at  F  to  the 
force  and  moment  at  S  for  the  steady-state  problem.  * 

Dynamic.  Now  consider  the  two  frames  in  their  •  lically  perturbed  (F"  and  S") 
states  (tig.  8).  The  perturbed  position  and  orientation  are  related  by 


RS"S'  +RS‘F'  +&F'F"  +RF"S"  =0 

(7. 1.1-8) 

where 

nF"  S"  T?F'  S' 

nS"  —  zts  i 

(7.1. 1-9) 

and 

CS"S'CS'F'CF'F"CF"S" 

(7.1.1-10) 

Thus, 

rf:f‘  =  cf"s"(rss:s'  +  rf::s")  -  cf"s"cs"s'rf;s' 

(7.1.111) 

pH  nl  r«/l  /'({  £»f#  C*  pl 

qF  r  _  qv  S  b  qS  F 

Taking  the  variation  of  both  sides  yields 

sRpr  =  cf"s,,(6rss:s'  -  6cs"s'rfs:s') 

c"C  -s"  S' 


■  S"  S' 


=  Cf"s"\6Rss:s'  +6*1."  (A -*5,  )iZ?;s'| 
St&r  -  Cf"s"64,§:s' 


(7.1.1-12) 


The  transformation  from  the  F"  frame  to  the  5"  frame  can  then  determined  in  terms 
of  R.  and  KG. 


% 


g F" S"  _  qF"S"  rF'S' 


(7.1.1-13) 


0  Cr”5' 

where  the  columns  of  Tt  are  associated  with  variations  of  the  generalized  coordinates 
JZf„s  and  8tl)g„s  ,  and  the  rows  are  associated  with  6Rp„F  and  6ipp„F  .  Then, 


Kl 


0 


0 


0  f£r£.s 


(7.1.1  14) 


where  the  columns  of  KG  are  associated  with  generalized  coordinates  Rf„s  and  6f,,  9  , 
and  the  rows  are  associated  with  6R^„S  and  6  V’,? „s  .  Equations  (7.1.1-13)  and  (7.1.1 
14)  define  the  constraint  formulation  for  dynamic  perturbations  about  the  steady  state 
configuration. 
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7.1.2.  Structural  Node  Demotion 


The  structural  node  demotion  constraint  describes  a  rigid  connection  1  tween  a  node 
D  in  the  child  subsystem  which  has  a  frame  of  reference  F ,  and  a  node  /  in  the  parent 
subsystem  which  has  a  frame  5  (fig.  9).  The  degrees  of  freedom  for  the  dependent  node 
D  will  be  eliminated,  while  the  degrees  of  freedom  for  the  independent  node  /,  the  super 
ordinate  frame  5,  and  the  current  frame  F  will  be  retained.  In  GRASP,  this  constraint  is 
generated  internally,  and  is  not  available  through  the  user  interface. 


i 


Figure  9.  Structural  node  demotion  constraint. 


Steady  State.  The  governing  equations  express  the  displacement  nnd  orientation  of 
node  D  in  terms  of  those  of  /,  5,  and  F.  The  basic  equations  for  the  deformations  come 
from 


R 


D'D 


R 


D‘V 


+  Rri  +  R 


/s'* 


rS’S  +  Rst  +  Rft  +  Rr 


FF' 


,F'D 


pn'l)  _pL)‘l’pri(plS,pS‘ScSFpFF'£F'D 


(7.1.2  1) 


In  matrix  notation  the  basic  equation  governing  displacement  is 

CFSCSS' \Cs'l(Cu' Rf,'1'  4  R1,'1)  t  R1/  4 


RD'D  pl>F'  fF’FnFSr’SS’  \  nS’lmll'  nD‘ /' 


Cnl  ’{Cri  Rfrs  *■  Rlf,1'  4  R$r ) 


F’F  ,  nOF’ 


(7.1.2  2) 
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It  is  also  necessary  to  take  the  force  and  moment  at  D1  and  find  their  contributions  to  the 
forces  and  moments  at  F\  and  S'.  For  this,  the  virtual  displacement  and  rotation  of 
D1  relative  to  D  is  required.  The  virtual  displacement  is 


,  ,~F'F 


6R°'d  -  -  Cdf'Wf,  CF'FCFSCss'[Cs,I(ClFRf;F  +RFT)  +  RIsf  +Rl',s}  + 


cue‘cF'rcFScss'HSs,s{cs,l(c,I'n?:,'  +  B.i;i)  +  Rf+R%s]+ 

CDF'CF'FCFSCSS'  [Cs‘\H'i rc"'RP'''  +  SR'/')  +  + 

CDF‘Si,l'FCF'FCFSRFS  -  CDF'6RFF 


=  -  Cdf‘6r£f  +  Cds'SRss:s  +  cd,sr’;'- 


cut  sty  r%f  4-  cDS  6ii>y  r$s  +  cuis^YRi 
~-cdf'(6rff:f  -Ag;F^f:r)+ 

Cds'(6Rss:s  -  &gsS*p)+ 

CDI(6RFl  -  Rf'^SrJ,1;1) 


D's  ,  nDI'i 


and  the  virtual  rotation  is 


(7. 1.2-3) 


6tj>s°  =  cdis^;1  +  c°*6+$s  -  cDF'srF:F 


(7. 1.2-4) 


The  virtual  work  performed  in  D'  contributes  the  following  terms  at  I',  S',  and  F' 

SW  »  -  (in£F)rF°:  +  (6R$,sf  Fg'  +  (SR‘i‘)‘Ff  -  (S^'Ff  (M$  +  R%fF°')+ 
+R$sFg,)  +  (H'I'ljr(Ml’'  +R?'FFl>') 

(7. 1.2-5) 

Dynamic.  For  dynamic  perturbations  about  the  steady-state  configuration,  the  basic 
kinematic  relation  is  used  to  determine  the  matrices  R  and  Ka .  The  basic  equation 
governing  the  displ*  ;emcnt  is 


BD"n>  —fiD”  l"  +  Jji-r  +  RVi  +  Ris"  +  RS"S>  +  RS'S . 
RFS  -  RF  F  -  Rf"f'  -  Rdf"  -  rD'd 


Rn"0'  zzBD'T  +  Sfl'  +  fi/'I  +  JJ/5"  +  ^5-5-  _  RF'S'_ 
rf"f>  rdf"  _RD'n 


(7. 1.2-6) 
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Expressing  the  position  of  the  perturbed  state  relative  to  the  steady-state  position  in  the 
D  basis  (fig.  9),  the  dependent  node  displacement  may  be  obtained. 


-,CDF  CFF  CF  S  Cs  5  Cs  lCu  Rfn1  + 

Cdf"cf"f'cf's'cs's"cs"i{ri"i'  +  Ri'ih 

CDF"CF"F>CF'S>CS'S"[RIS"  +  Jtsy'y 

qDF"  (jF"F'  qF'S'  RFj  S' 

rDF"  /  dF"  F1  ,  dDF"  \  d D'D 
v  \Jn-Fn  T  -tf-pn  )  — 


(7. 1.2-7) 


The  first  variation  is 


-CDF"  CF"  F>  CF>  s>  Cs>  s"  Cs" 1 8R*"  -f 
cdf"cf"f'cf's>cs>s"srs';s'  _  cdf"6Rff;f'  + 

cdf"cf"f'  cf's'  cs's"  Cs" 1 C11’  8C1'  iU  R^" 1,1 4- 
CDF"CF"F>CF>s'scs's"[cs"rcn>ci'i"Rv;'"+ 

C*'I(Ri't'  +R1I'I)+  R^"+Rss:s'}  + 

CDF"6CF"F'CF'S'  [cs's"cs"icii'ci'i"rd';i"+ 

CS'S"CS"I{RI"I'  +  RVI)  +  cs's"(4f  +  Rss:s')  -  Rf:s '} 


(7. 1.2-8) 


Similarly,  the  relationship  governing  the  virtual  rotation  is 

W_D"D'  =^”/"  +  6J,1"1'  +6J,1'1  +  6JLIS"  +  8JS"S'- 
SjS's'  -61f"f'  -8J,df"  -8±d'd 

which  in  the  D  basis  becomes 

6^"°'  =CDF"CF"F'cF's'c5's''C5''/^f/'  + 

cdf"cf"f'cf’s'cs's"hs:s'_ 

Cdf"8Vf:f' 


(7. 1.2-9) 


(7.1.2-10) 


39 


The  virtual  work  done  by  a  force  at  D"  is  then 


(6R%'D')TFg"  =(6Rf‘')T(F,D'  +  c,*pge£3'  -  C"”f£'<Cf’)+ 

tc  TiS"  S’  \T (  J?D’  ,  £,D'  A S"  S'  nS’F'  fi D'  oF"  F'  \ 

{OJtgn  )  [±$1  *r  ^5/  Ugn  —  O  2*jri  i/jp/i 

(«#:'”  )T*#+ 

W,,",')T(R?','F?'  +  Ff'Rf'Fei"1'  +  c,s'R$s'PB'e  p'- 

cif'r$‘'fi?:6f:f')+ 

{H>£srfiRg*Fg  -  CS',FP'R‘I",‘  -  F0RsP  + 

CF‘Fl>'Rf>' +  Rp'F fep'- 

rs’F 1  r>D's'  irD'aF"F'\ 

O  Jttpi  r  pi  vpn  )  — 

(h£:f')t(.r?;f'f!?:  -  cf'' fp'r p  -  cF's'FpRp'+ 
cF',Fp'R?',,sf’‘ + cf's'f§‘ Rg:s'e§::s'- 

nD1  f)D' F'  nF"  F1  \ 

Jr  fi  it  pi  vp'"  ) 


and  the  virtual  work  done  by  a  moment  at  D"  is 

wg"D,)TMg"  =wp'fc,s"(  a  -  si  :s')cs'f\& + e£,F')M?::+ 
(«*!,?■  f(A  -  ei :$')cs'f'(a + eF::F’)M°::- 
wp'fMp: 

=W-,"'')T(MP'  +  cIS'M°'e i:s'  -  c'F'MpeF';F')+ 


(6i,p')r(Mp  +  Mpep'  -  cs'rM$e p~) 

(SVP'fMp, 


,S'F'  fjD'  nF"  F' ' 


(7.1.2-11) 


(7.1.2-12) 


The  matrices  R  and  Kc  are 


c°i  -CDIRf‘F  CDS'  -Cds'R$s'  -Cdf'  Cdf'R$f' 


CDS' 


- CDF ' 


(7.1.2-13) 


where  the  columns  of  %  are  associated  with  variations  of  the  generalized  coordinates 
SRj",1' ,  Sijif  ',',1' ,  8R^"S',  8‘ifsg", s' ,  6Rp,',F' ,  and  8il)F",F' ,  and  the  rows  arc  associated  with 
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SRp",0'  an dSipp,',0'.  Then, 


K 


G  _ 


0 
0 
0 


0 

f  -  rV 


-fP'r?'1' 


Cs'1Ff>'  -C^tFp'R?'1' 


-CfiFP'  CF1 FP'  Rf'1' 


0 

0 

0 

pp: 

0 

-cf's'fp: 


o  ciF'pp; 

o  cif\r»:i'fp: +m$) 

0  CSF'  Fp! 

0  C** {R$s' FP'  +  Mg) 
0  0 
0 


where  the  columns  of  K G  are  associated  with  perturbations  of  the  generalized  coordinates 
Rj",1' ,  Oj',',1',  Rs'n5' ,  Os”' S' )  Rp'i'iF' ,  and  0p!F',  and  the  rows  are  associated  with  SRp,1' , 

w:,1',  6Rss:s',  s^:s\  srf:f',  and  6^;f\ 


7,1.3.  Screw 

The  screw  constraint  describes  two  nodes,  D  and  I,  that  are  connected  by  a  mecha¬ 
nism  that  permits  translation  along  and  rotation  about  a  single  axis  which  is  fixed  in  the 
coordinate  system  of  both  nodes.  The  dependent  node  D  will  have  its  degrees  of  freedom 
eliminated,  while  the  degrees  of  freedom  of  the  independent  node  I  wil'  be  retained.  This 
constraint  is  available  through  the  GRASP  user  interface. 

To  simplify  the  derivation,  two  intermediate  nodes  located  on  the  screw  axis  e>cr  will 
be  introduced  (fig.  10).  The  “stationary”  node  5  is  rigidly  connected  to  the  independent 
node  I ,  while  the  “moving”  node  M  is  rigidly  connected  to  the  dependent  node  D.  Nodes 
ill  and  5  initially  coincide  in  both  position  and  orientation,  but  may  translate  along  and 
rotate  about  the  screw  axis  relative  to  one  another. 
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These  equations  simplify  somewhat  since  CIS  =  Cl  S  —  CD  M  =  CDM  =  A.  CM 
can  be  easily  expressed  as  an  Euler  rotation,  given  the  screw  rotation  8i  and  effi  ,  the  screw 
axis  unit  vector.  The  virtual  displacement  is  then 

5R$d  =CDISCn'  Cvs‘  C*VM'  CM'D'  i2g:w'  + 

qDIqII'  gf'S'QS'lU'  CM'D'  R&m'  + 

CDI8CTI>  Cls'  eg/  u  +  CDICII'cI's'c‘/Su+ 

ca,(sc"‘Rf.''  +  sn'i')  t7-1-3'4) 

=CDI  [U if'  -  ( Rf‘ +  u er'  +  Rf 

ctIcr'8u  +  e'Jcr'R?'M'60\ 
and  the  virtual  rotation  is  given  by 

6$%°  =  CDI(8ij)jI  +  e'Icr'  80)  (7. 1.3-5) 


The  virtual  work  at  the  screw  connection  and  at  /'  due  to  a  force  and  moment  at  D' 


is 


SW  =(8Rjr)TFP’  4-  (Rf  1  +  Mf  )+ 


i'i 


T.  - 


D'  I'  i pD' 


rD'< 


Me/crV  i7/3'  +  60(etIcr'  )T  (Rf  M>  Fp'  +  M f) 


(7.1.3-6) 


Dynamic.  For  the  dynamics  of  the  constraint,  the  equations  governing  the  degrees 
of  freedom  arc  used  to  to  find  the  matrices  7 Z  and  K G.  Consider  the  nodes  and  the  screw 
axis  in  their  perturbed  states  (fig.  10),  an  infinitesimal  perturbation  from  their  steady-state 
positions  and  orientations.  The  basic  equations  are  similar  to  those  of  the  steady-state  case. 


=:RD"M"  +RM"*<"  +Rlit"S" 

ri",,  +  ri'd' 


+  RS"J"  + 


6jpD"D'  =6jD"M' 

Sil*5"1"  +8rj,J"s'  +8il>1,u' 


-I-  8tp 


M"Sl" 


+  8i ’> 


4- 


(7. 1.3-7) 


The  first,  third,  fourth,  and  sixth  terms  arc  zero  in  both  equations.  Proceeding  as  above, 
and  noting  that 


Cl  1  =A  -  Cri0II"I'C!l'  =  CI!(&  - 

i"i' 


8Cr  1  =  -  Cr,( A  -  0rr)8ii  '  C 


,i’i 


,W 


8C""m"  -  b0?$,'CM"ht" 


(7. 1.3-8) 
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the  7t  and  KG  matrices  are  then 


77  = 


CDI  -CDIkyr  CDIe'fr'  CDle?r'R?'M' 
0  CDI  0  CDIetjcr' 


(T.l.3-9) 


where  the  columns  of  7 Z  are  associated  with  the  variations  of  the  generalized  coordinates 
6Rj  1  ,  Srpj  1  ,  6u,  and  66 ,  while  the  rows  are  associated  with  8Rq  d  and  8tJ)q  d  . 


KG  = 


0  0  0  0 
0  -Pp'R?'1'  -ej-'if'  -Ff'Rf'M'er' 


0  -(e/cr,)T-P/y 


0 


0 


(R?'w')T«r'Ff'  (RD'M1)1  W  per'  pD' 

°  flit?'  °  -(■«/)  eI  eI  *1 


(7.1.3-10) 


where  the  columns  of  KG  are  associated  with  the  perturbations  of  the  generalized  coor¬ 
dinates  R\  1  ,  6\  1  ,  u,  and  6  and  the  rows  are  associated  with  6Rj  1  ,  6ipj  1  ,  6u ,  and 
66. 


7.1.4.  Copy 

The  copy  constraint  describes  the  relationship  between  generalized  coordinates  that 
are  common  to  both  parent  and  child  subsystems,  but  are  otherwise  unconstrained.  This 
situation  most  often  exists  when  unconstrained  generalized  coordinates  in  the  child  sub¬ 
system  are  passed  up  to  the  parent  subsystem.  This  constraint  is  not  available  through 
the  GRASP  user  interface. 

Steady-State.  From  equation  (7.1-1),  the  constraint  relationship  between  the  child 
subsystem  generalized  coordinates  qei  that  will  be  eliminated  and  the  parent  subsystem 
generalized  coordinates  qn  that  will  be  retained  can  be  written  as 

,  (i  =  l,...,Ar)  (7 .1.4—1 ) 

Therefore,  the  calculation  of  the  contributions  of  these  generalized  coordinates  to  the  child 
subsystem  state  vector  involves  only  copying  the  values  of  generalized  coordinates  from 
the.  parent  subsystem  state  vector  into  the  child  subsystem  state  vector. 

The  variation  of  qe.  is  then 

(i  =  1,. . .  ,  AT)  (7. 1.4-2) 

When  this  expression  for  6qei  is  substituted  into  equation  (7.1-2),  the  virtual  work  is 

N 

<W  =  £*3,.,  (C?r, +(?«,)  (T.l.4-3) 
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Thus,  to  assemble  the  residual  vector  for  the  parent  subsystem,  the  contributions  of  these 
generalized  forces  from  the  child  subsystem  are  added  to  the  generalized  forces  from  the 
parent  system. 

Dynamic.  The  derivation  of  the  constraint  dynamics  follows  a  similar  vein.  First, 
the  perturbed  generalized  coordinates,  the  variations  of  the  generalized  coordinates,  and 
the  generalized  forces  are  expanded  in  Taylor’s  series. 


—  Qrj  j 

% 

t-T 

II 

■<s> 

(7, 1.4-4) 

8qCi  =  8qrj 

(i-l . JV) 

(7. 1.4-5) 

N  N 

Qe;  =Qei  +  ie"V  *7  » 

j=l  ;=i 

N  N 

Qri  —Qn  4  ^  ]  Lriej  qCj  +  Lr;rs  qr. , 

i=  i  i=  i 

When  these  expressions  are  substituted  into  e(piation  (7.1-2),  and  the  resulting  expression 
simplified,  the  virtual  work  is  written  as 


(7.1.4-6) 


(i  =  1,...,JV) 


N 

8W=Y,6qri 

i-l 

Since  the  qe.  generalized  coordinates  exist  only  in  the  child  system,  and  the  qri  exist  only  in 
the  parent  system,  Le(ri  and  ~L*.r.  are  null.  The  7 Z  matrix  is,  therefore,  an  identity  matrix. 
For  small  perturbations  about  the  steady-state  solution,  the  coefficients  in  the  rows  and 
columns  associated  with  the  copied  generalized  coordinates  in  the  child  subsystem  dynamic 
matrices  (Hi,  C,  and  K)  are  simply  added  to  coefficients  in  the  corresponding  rows  and 
columns  of  the  parent  subsystem  dynamic  matrices.  The  geometric  stiffness  matrix  KG  is 
null. 


N 


Q,:  +  Q«  +  £(£  ctej  4"  f'ejr/  4*  Lrie-  +  Lr{r. )  qr. 
7=1 


(7. 1.4-7) 


7.1.5.  Prescribed 

The  prescribed  constraint  is  used  to  describe  the  permanent  deformation  of  a  partic¬ 
ular  generalized  coordinate.  This  constraint  is  trivial,  because  the  steady-state  value  is 
constant.  In  GRASP,  the  prescribed  constraint  is  available  through  the  user  interface  for 
nodal  degrees  of  freedom. 

Steady-State.  Following  the  derivation  of  a  general  constraint,  consider  a  child 
subsystem  that  has  Nc  generalized  coordinates  qe.,i  -  1,...  ,  iVe.  For  this  constraint,  one 
of  those  generalized  coordinates  ( e.g qCx)  has  a  prescribed,  constant  value. 


qei  =  constant 


•7c,  —  <7:  ((Jr,  i  ■  •  •  >  QrNr  ) )  (i  —  2, ...,  N  ) 


(7. 1.5-1) 
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The  total  virtual  work  is 


N' 

SW=Y,heiQ', 

i-i 
Nc 

=  J26qeiQci 

i—2 

since  the  variation  of  the  prescribed  generalized  coordinate  Sqe i  is  zero.  Therefore,  this 
generalized  coordinate  makes  no  contribution  to  the  virtual  work  of  either  the  child  or  par¬ 
ent  subsystem.  In  practice,  degrees  of  freedom  that  are  prescribed  in  the  child  subsystem 
may  be  eliminated  from  the  parent  subsystem  state  vector. 


NT 

53  SqriQn 

i=i 

Nr 

53  V»Qr, 
«=1 


(7. 1.5-2) 


Dynamic.  The  derivation  of  the  dynamic  constraint  equations  for  small  perturba¬ 
tions  about  the  steady-state  solution  proceeds  following  equations  (7.1-5)  thiough  (7.1-7). 
The  only  difference  is  that  in  equations  (7.1-5)  and  (7.1-6),  i  —  2, ...,iVe.  When  these 
expressions  are  substituted  into  equation  (7.1-2), 


Ne  AT 

<w=EXX 

t=2  j=l 


dgi 

dqr. 


d2gi 


NT 

+  S  0q 


ru 


Nt  Nr 

■4rh  +  *  •  •  j  1  Qa  +  53  ^e<e»,9efc  H*  E  ■^ejrfc  Qri, 

k-1 


k= 2 


Nr  /  N'  Nr 

+  53  Sqr<  I  +  53  Lr>e>q*i  +  LuriQr; 

i=  1  \  j—2  j=  1 

(7. 1.5-3) 

From  this  equation  it  can  be  seen  that  the  contributions  to  the  virtual  work  are  the 
same  as  for  the  general  case,  with  one  exception.  The  rows  of  Le.eh  and  ie;rfc ,  and  the 
columns  of  Leieh  and  Lr.Cj,  associated  with  the  prescribed  generalized  coordinate  have 
been  eliminated.  This  is  equivalent  to  removing  the  appropriate  rows  and  columns  from 
the  M,  C ,  and  K  matrices  that  are  passed  up  to  the  parent  subsystem. 


7.1.6.  Copy  Air  Mass 

The  copy  air  mass  constraint  is  the  constraint  used  to  transform  the  air  mass  gener¬ 
alized  coordinates  and  forces  between  cltild  and  parent  subsystems.  This  constraint  is  a 
clone  of  the  copy  constraint,  specialized  to  copy  only  the  four  air  mass  degrees  of  freedom. 
Due  to  the  fact  that  the  air  mass  degrees  of  freedom  aie  defined  in  an  inertial  frame,  and 
need  never  be  tranoiormed  out  of  that  frame,  the  generalized  coordinates  and  forces  are 
simply  copied.  The  copy  air  mass  constraint  is  not  available  through  the  GRASP  user 
interface. 


7.1.7.  Periodic  Frame 

The  periodic  frame  constraint  describes  the  relationship  between  a  superordinate  (par 
ent)  frame  5  and  three  or  more  identical,  child  frames  Ff.  (for  k  ~  1,2,..., 6)  rigidly  at¬ 
tached  to  5.  Frames  Fk  are  located  at  equally-spaced,  azimuthal  intervals  about  an  axis 
fixed  in  S  (fig.  11).  The  origin  of  5  is  located  on  the  axis  of  symmetry,  while  the  origins 
of  the  Fk  may  be  located  elsewhere.  In  GRASP,  the  periodic  frame  constraint  is  available 
through  the  user  interface. 
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Figure  11.  Periodic  frame  constraint. 


The  derivation  of  the  periodic  frame  constraint  is  very  similar  to  that  for  the  fixed 
frame,  except  that  it  is  assumed  here  that  there  are  b  identical  frames  spaced  at  equal 
azimuthal  intervals  around  an  axis.  The  quantity  Rp?k  is  independent  of  k,  and  CFkS  - 
CFlSTk  where 


Tfc  =  T0  +  Tc  cos  (pk  4-  T,  sin  (f>k 


(7. 1.7-1) 


and  where 


To 


10  0 
0  0  0 
0  0  0 


Tc  = 


0  0 
0  1 
0  0 


0 

0 

1 


Ts 


0  0  0 

0  0  1 

0-1  0 


<f>k  =  -y(fc  -  1))  fc  =  l,2>...,6 


(7. 1.7-2) 


(7. 1.7-3) 


The  fixed  frame  equations  ran  he  easily  modified  to  account,  for  this  configuration.  It  will 
he  assumed  that  the  axis  of  symmetry  for  the  periodic,  child  frames  is  . 
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Steady- State.  For  the  steady-state  problem,  the  equations  for  the  deformed  position 
and  orientation  of  any  one  of  the  child  frames  can  be  written  as 

RFj,Fh  -~-.CFis‘TkRs's  +  RFj,S'  -  CF'»FkRFkS 


=CF's'TkRss',s  +  RFJlS‘  -  CF^s,TkCs'sT^CSFlR^s 


(7. 1.7-4) 


CF'kFk  =CF'>s‘TkCs'sT?CSF' 

In  order  to  make  the  left-hand  sides  independent  of  fc,  let  i?J,f  =  =  0  and  (/>f  5  = 

(j) f  6  —  0.  Since  the  right-hand  sides  are  equivalent  for  all  k,  all  Tk  can  be  set  to  T0  to 
simplify  the  equations.  The  virtual  displacements  are  then 


6RFj,Fk  =C^s'To6R%s  -  CF‘^S'T()6CS'STFCSF^R^S 


=Cf's't08Rss‘,s  +  CF'S' Tq6iI>Ss,S Cs  sTn  CSFl RpkS 


(7. 1.7-5) 


and  the  virtual  rotations  arc 


tyS*  =  cF^s'wi',s 


(7. 1.7-6) 


The  virtual  work  at  S'  due  to  the  b  sets  of  forces  and  moments  acting  at  F l  is  therefore 


6W  =  j2[(6RnFkf  FfI  +{WFJlFk)TMFl 

k=  1  L 

~b^(6R%s)T TFCS' F'1  F Fj  4-  (^f!5)7  [TFCs'F''MFi 

{tfcs'f'  FFi)~cs'sTQrcSFi  f?£s] } 


(7. 1.7-7) 


Dynamic.  For  small  perturbations  about  the  steady-state  solution,  the  perturbed 
position  and  orientation  (fig.  11)  of  any  one  of  the  child  frames  Fk  is 


RFFk„Fk  +CF"s"TkRSs,',s'  -  CF" s" TkCs" s' Tk  Cs' F'1  Rpk,S 


Fit  iTi I  n't*  C’ '  c't  o*  m  c*  pi 

hph  -cF* 5  Tkcs  stTcsf 1 


(7. 1.7-8) 
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To  first  order  in  the  perturbation  quantities,  the  virtual  displacements  and  rotations  are 


61#*  =CF"s"Tk  +  ^s7( A  -  fs„S  )T?CS'*R*S' 

rh  * 

Fn  I?*  t e"  c*1  c1 

6 A*  =C*S  Tk6il>§„s 


(7.1.7-9) 


-s  s 


Note  that  the  geometric  stiffness  matrix  will  come  from  the  9S„  term. 
After  substituting  into  the  expression  for  the  virtual  work,  the  matrix  77  is 


71  = 


CF"s"Tk 

0 


-(T’FCs'FiRnS'Fiy 

CF"s"Tk 


(7.1.7-10) 


where  the  columns  of  77  correspond  to  6Rf,',s'  and  firf)§"s' ,  and  the  rows  correspond  to 
6R*„*  and  bill*,* .  The  matrix  KG  is 


Kg  = 


0  0 
0  j2T:C*KF*R*S'cWTk 


k=  1 


(7.1.7-11) 


where  the  columns  of  KG  correspond  to  Rg„s  and  0f„5  ,  and  the  rows  correspond  to 

pill  pit  pill  pil 

6RJi,  h  and  fiipr/,,  ' .  For  evaluation  of  the  lower-right  submatrix,  it  should  be  noted  that 
r  ** 


b 

Y  Tk(  )H  =  H r,T(  )H  +  }2?(  )%  +  \Tj(  )T,\  (7.1.7-12) 

k-\ 

when  the  expression  enclosed  in  parentheses  is  independent  of  k. 

7.1.8.  Periodic  Node  Demotion 

Just  as  the  periodic  frame  constraint  is  very  similar  in  concept  to  the  fn:ed  frame  con 
struint,  the  periodic  node  demotion  constraint  has  a  similar  relationship  to  the  structural 
node  demotion  constraint.  In  this  case,  a  node  belonging  to  a  parent  subsystem  is  repli 
cated  in  the  child  subsystem  at  b  equally  spaced  azimuthal  intervals  about  an  axis  that 
is  fixed  in  the  parent  subsystem.  The  periodic  node  demotion  constraint  is  not  available 
through  the  GRASP  user  interface. 

The  degrees  of  freedom  of  the  b  child  subsystem  nodes  Dk  are  expressed  in  terms  of 
the  degrees  of  freedom  of  Fk,  S,  and  I.  To  visualize  this  constraint,  consider  figure  11  and 
imagine  a  node  I  associated  with  frame  5  and  a  node  Dk  associated  with  each  frame  Fk, 
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as  in  figure  9.  The  virtual  work  done  at  all  nodes  Dk  for  k  =  1,2, . . .  ,b  is  determined  at 
Fk .  The  total  virtual  work  is  summed  for  k  =  1,2, . . .  ,6  and  determined  at  S  and  I.  For 

this  constraint,  it  is  assumed  that  the  axis  of  symmetry  is  ,  and  that  also  lies  along 
that  axis. 

Steady- State.  The  governing  equations  for  the  periodic  node  constraint  are  derived 
in  a  similar  manner  to  those  of  the  structural  node  demotion  constraint.  First,  let 


CF'S  =CF*S'  =  CFlSTk  =  CFis'Tk 

CDhI  =cO!  /'  _  CD,ITk  =  CD[I‘Tk 
where  Tk  is  defined  in  equation  (7. 1.7-1)  and  (7. 1.7-2).  When 

—Rpk,S  =  constant 

R^J  =R^kI  —  constant 

the  positions  and  orientations  of  the  frames  and  nodes  may  be  written  as 

R^Dk  =CDhF'kCF'kFhCFlSTkC SS'  + 

izfr)  +  -R^'  +  ^5]- 

C^n(CF^RFf  +  R$f'R%?Fk) 

CD'"Dl'  =CD'I'TkCI'ICIS'cs'sT?CSF'CFkF'''CF!'Dk 


(7. 1.8-1) 


(7. 1.8-2) 


(7. 1.8-3) 


To  make  equation  (7. 1.8-3)  independent  of  k,  let  Rg?  =  0,  CIS  =  A  and  let  Tk  =  To  (all 
displacements  and  rotations  for  S  and  I  take  place  along  or  about  the  axis  of  symmetry). 

The  virtual  displacements  and  rotations  are  required  in  order  to  calculate  the  virtual 
work  of  a  force  and  moment  at  D'k  for  all  k. 


SR'& 


n'kD’'  =  -  CD''F'1' 6i/>Fpfh CF,iFk CFlSTvCss'  \CS>  1  {C11' T$ CIDl 7 ’  + 

R\F^Risr +rss:s]+ 

CDkF'k CF>hFk CFiSTqCss' Srj}Sg,S  [cs> i (Cn' Tq  CIDl Rp’!'1'  + 

R'ft+R*/ +rss:s}+ 

CDhF'hCFlFhCFxSncSS'  {cS'l(^1ClFTFCI,D'>R%f  + 

SRFl)i-6Rss:s}- 


CDkF'>{-  b/F\ F"  C F ^ Fk  RFl s  -l-  6 RFj, F" ) 


(7. 1.8-4) 


Fk 
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(7.1.8-5) 


=c°'FLcFlF'cF'sT<,css'{cs',stf'‘  +  S4$s)~ 

CD'F*WFlF' 

The  virtual  work  due  to  the  virtual  displacements  of  each  of  the  nodes  D'k  is 


«W  =  -  (6RFJlFk)T F*  -  (6tlipj'Fk)T(Mp'?  +  R*FkF*)+ 
(6Rss',sfcs'sT?CSF'Fgk+ 

( Ws'  S)T{CS'  sT0r  CSFl  Mpk  +  [Cs'sTFCSFlCFkF!‘(Rp}Fk  +  CF^FkRtfks+ 

CFi,F>CF'sT0Css,Rss:s)  ]  "Cs'sTFCSF'Fp}  }  + 

(8Ri;i)TCIS,Cs'sTFCSF'F?y\- 

(, Srj>Ii1)T{CIS'cs'sT^CSFiMpk  +  [CIS‘ Cs  sTFCSFl CF|,F*  + 

CnF»RF>s  +  CF*FkCF'sT0Css'{Rps  + 

CK  fh  cfi  sToCss'  (js'i^'i} }  “  cis'  cs'stt  csFi  Fpl ) 

(7.1.8-G) 

The  summation  cf  terms  involving  the  virtual  displacements  and  rotations  at  S'  and  /' 
involves  only  a  multiplication  by  6.  The  corresponding  terms  at  Fk  need  not  be  summed 
since  only  one  system  contributes  to  the  virtual  work  there. 


Dynamic.  As  in  the  case  of  the  equations  for  the  steady-state  periodic  node  demotion 
constraint,  the  derivation  of  the  equations  for  the  dynamics  is  similar  to  the  derivation  of 
the  structural  node  demotion  constraint  dynamic  equations.  In  a  manner  similar  to  that 
for  the  static  equations  above,  let  Rtf  =  0  and  CIS  =  A.  Also,  let 


CFhS  =CFlS'  =  CFl'S"  =  CF,STfc  =  CF[S'Tk  =  CF^"Tk  =  constant 
CDkI  =CD'<'i"  -  CD^"  =  CDlITk  =  CD'^Tk  =  CD"rTk  =  constant 
RFtf  =RFkS'  =  RF!is"  =  constant 
iZy* r  =Fy|' ;  =  R^y,1  =  constant 


(7. 1.8-7) 
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The  resulting  equation  for  the  matrix  71  is 


n  = 


CDkI  -R%kI'cDkI  C DkS’  -CDkS>  R^!'S' 


0 


CDkI 


0 


CDkS' 


_cd>k  cDkFiR^F 


n 


- c °>>n 


(7. 1.8-8) 


where 


CD„I  =cD-nCF[s'TkCs'r 


CDhS‘  =CDkFkCF'^s'Tk 

R*s‘  =Cs'Fh?CFlD‘R%''  +  Cs''Rfi  +  R'/ 


(7.1.8-9) 


The  columns  of  71  correspond  to  SRj  1  ,  6ipj  1  ,  SRss,  5  ,  6tp^ns  ,  6Rp„  h ,  and  k, 


\S" s'  tnn'n  i  cFL'fl 

*  t  t  D"  D1  Du  D1 

respectively;  while  the  rows  are  associated  with  h  and  6if>D*  h.  The  coefficients  of 

the  geometric  stiffness  matrix  KG  are  then 


5iJj  1  row: 


r\In  I1  i 

it}  column:  0 


Of  1  column:  0 


Rs"S  column:  0 


0'l„s  column:  —b{TjFfx) 


RhkF"  column:  0 


(7.1.8-10a) 


OpnFk  column:  T^CIFlFp,x 


f: 
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C  ,  I"  I 

o-ipj  row: 


Rj  1  column:  0 


0}"r  column:  -  bT*{F?1  R?lI)T0  -  J-Tj{F^  R^ll)l 

M 

~Tj(F?'R°'',)T. 


b-T,T?D[6D[I, 


Rs',S  column:  0 


6s5"s'  column:  -  bT^R?'1' F?'l)T0  -  -T?(R?[I  F?x 

h-Tj(R?','F?')T,  -  b(TjM?'')~ 


Rph„Fk  column:  0 


eFlK  column:  T?  R?[I>  CIF'Fp!1  +  TkCIF'M$ 

rh  1  rl  rl 


5Rt 


S" 


row: 


Rj  1  column:  0 


6{  1  column:  0 


ilfii5  column:  0 


#f„5  column:  —  ^T^F^1) 


Rpn Fk  column:  0 


dp1;, Fh  column:  TkCIF'Fp,1 


(7.1.8-10b) 

)Tc- 


(7.1.8-10c) 
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6ipgnS  row; 


Rj  1  column:  b(T{f  Ff*1) 


column:  -  M? '(#f*  k“'‘  )T„  -  ^TcT(F,D‘RfJ')T<- 
~T?(FID''a?‘,')T, 


Rs,',s'  column:  b{Tj F^1) 


ess:s'  column:  -  F?')T0  -  ”-T?(R?'*  F?')TC- 


b~T,»D[l'  i~D\ 


Ffl)T.  -  b(T0TMp) 


rt II  r<< 

R'Jl  fc  column:  0 


0% K  column:  Tfcr R?[1' CIF>' F°!'  +  TkCIF' 


rptl  pf 

6Rp„  h  row: 

Rj"1  column:  0 
Oj  1  column:  0 
Rg„ s  column:  0 
8g„s  column:  0 
Rpi,Fu  column:  0 

9~nFk  column:  0 


(7.1.8-10d) 


(7.1.8-10e) 
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8ij) 


i,K'K 


K' 


row: 


Rj  '1'  column:  -  F^CF'^Tk 
e'j"1'  column:  FpfC^1  R?'1' Tk 
Rss",s'  column:  -  F$CF'ITk 
ef,',s'  column:  FF!lCF'^Rf'lI'Tk 


(7.1.8-10f) 


pit  pi 

Rpn  k  column:  0 

aFi'FL  i  BDinniFi 

9 pi,  *  column:  -  Ff,  Rf, 

7.1.9.  Periodic  Generalized  Coordinate 

A  rotationally  isotropic  structure  consists  of  three  or  more  identical  substructures 
that  are  spaced  around  an  axis  of  symmetry  at  equal  azimuthal  intervals.  The  periodic 
generalized  coordinate  constraint  exists  in  order  to  transform  generalized  coordinates  that 
belong  to  the  rotationally  isotropic  structure  into  the  generalized  coordinates  for  a  generic 
member  of  that  stiucturc.  Additionally,  it  must  transform  generalized  forces  for  a  generic 
substructure  into  the  generalized  forces  for  the  complete  structure.  In  one  sense,  it  is 
simply  an  extension  of  the  copy  constraint  for  periodic  structures.  In  GRASP,  the  periodic 
generalized  coordinate  constraint  is  not  available  through  the  user  interface. 

Steady-State.  The  set  of  independent  generalized  coordinates  for  a  rotationally 
isotropic  structure  may  be  grouped  as  collective  (go),  cosine  (gc),  and  sine  (g,)  components. 
The  generalized  coordinates  for  the  kill  generic  substructure  qk  may  be  written  as 

Qk  =  go  H-  qc.  cos  <f>k  +  g,  sin  (f>k  (7.1.9-1) 

where  (f>k  =  —  1).  The  variation  of  these  coordinates  is 

f>qk  =  <5?o  +  6qc  cos  <f)k  +  8q,  sui<f>k  (7.1. 9-2) 

Given  generalized  forces  Qk ,  the  total  virtual  work  from  all  of  the  generic  substructures 

b  b  b  b 

<5W  =  V  8qkrQk  =  tqo1  Qfc  +  he1  X  ®k  cos<^  +  Sq>1  X  ®k  snx^k  (7.1.9-3) 

fc-J  k= 1  k=  1  k= 1 


IS 


Since  the  generalized  force  of  a  generic  structure  is  independent  of  k ,  Q0  =  Qk  and 

SW  =  6q0TbQ0  (7. 1.9-4) 

Dynamic.  The  dynamic  perturbations  of  the  generalized  coordinates  are  related  in 
the  same  manner  as  the  variations  of  the  steady-state  generalized  coordinates. 

qk  =  qo  +  qc  cos  fa  +  q,  sin  fa  (7.1.9-5) 

Like  the  generalized  forces,  the  substructure  coefficient  matrices  M,  C,  and  K  are  inde¬ 
pendent  of  k  and 

6W  =  6qJ{Mjk  +  c'qk  +  Kqk)  (7.1.9-6) 

The  contribution  to  the  virtual  work  in  terms  of  the  independent  generalized  coordinates 
is 


h 

6W  =  ^  [Sg0T  +  SqcT  cos  <f>k  +  ^g^sin 
fc=i  L 


fa\  |  (Mj 


q0 

qc  cos  fa 
sill  fa 


\  .  )  f  7°  ) 

[C\  <  gccos0fc  >  +  [K]  <  qc  cos (f)k  >  > 
{  qt  sin  cf)k  J  [  q,  sin  fa  J  ^ 


(7. 1.9-7) 


bdqo1  ( Mq0  +  Cq0  +  if  go)  +  -^gc^  (-^e  +  Cgc  +  Kqc) 

+  \h»T{Mq»  +  Cq.  +  Kq>) 

The  matrices  for  the  rotationally  isotropic  structure  therefore  have  three  rows  and  columns 
for  every  row  and  column  in  the  generic  substructure,  and  are  of  block  diago  nal  structure. 


7.1.10.  Periodic  Air  Mass 

The  periodic  air  node  constraint  describes  the  transformation  of  the  air  node  gener 
alized  coordinates  and  forces  between  subsystems  associated  with  periodic  structures  and 
subsystems  associated  with  generic  substructures.  Since  the  air  node  generalized  coor 
dinates  describe  an  induced  airflow  velocity  field  that  is  already  axially  symmetric,  the 
periodicity  of  the  structure  has  no  effect  on  them.  In  fact,  it  is  assumed  a  priori  that  the 
flow  field  is  interacting  with  a  rotating,  periodic  structure.  This  constraint  is  not  available 
from  the  GRASP  user  interface. 

Steady-State.  When  a  subsystem  is  periodic  (in  the  sense  that  it  consists  of  three  or 
more  generic,  periodic  members  such  as  those  described  under  the  periodic  node  demotion 
constraint),  the  steady-state  air  node  generalized  coordinates  Lr,4  and  ^/lr  are  simply  copied 
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from  the  parent  subsystem  to  the  child  subsystem  in  a  manner  similar  to  the  copy  air 
mass  constraint.  During  the  assembly  of  the  generalized  forces,  the  air  node  generalized 
forces  from  a  single,  generic  substructure  are  simply  multiplied  by  b  and  added  to  the 
corresponding  generalized  coordinates  of  the  parent  subsystem. 

Dynamic.  For  perturbed  motions,  let  the  generalized  coordinates  for  the  fcth  subsys¬ 
tem  be 

<ik  =  Tkq\  6qk  =  TkSq  (7.1.10-1) 

where  Tk  is  as  given  in  equation  (7.1. 7-1)  and  where 


and 


fii  ] 

r  PiA] 

1  ^i2fc  7  ; 

i 

q  =  < 

</4 

Itfl  j 

f  SP&  } 

1  i  SP,A  1 

6qk  =  < 

w* 

‘-o 

II 

{ J 

l  S+A ) 

(7.1.10-2) 


(7.1.10-3) 


Note  that  qk  and  q  will  not  appear  in  the  dynamical  equations.  The  equations  now  trans¬ 
form  in  exactly  the  same  manner  as  the  ones  in  the  copy  air  mass  constraint. 


7.1.11.  Rotating  Frame 

The  rotating  frame  constraint  describes  a  constraint  that  is  very  much  similar  to  the 

fixed  frame  constraint,  except  that  frame  F  is  rotating  at  a  constant  angular  speed  relative 

*  F 

to  frame  S  (fig.  12).  The  axis  of  rotation  passes  through  the  origin  of  F  and  along 
.  No  time-dependent  terms  are  retained  in  the  equations.  In  GRASP,  this  constraint  is 
available  through  the  user  interface. 

Steady-State.  In  moving  to  its  steady-state,  equilibrium  position,  the  axis  of  rotation 
*  F  _  C  P  E1* 

follows  6,  .  The  position  vectors  R  and  R*  are  constant  in  the  5  and  S'  bases, 
respectively.  The  change  in  orientation  is  then 

CF's\t )  =  T{t)CF's'{0)  =  CFS{t)  =  T(t)CFS{  0)  (7.1.11—1) 

where 

T{i)  =•■  To  +  Tc  cos (TP)  4-  T,  sin  (fit)  (7.1.11-2) 

where  I'o,  Tc,  and  T*  are  given  in  equations  (7. 1.7-2). 

Tlie  kinematics  for  the  rotating  frame  constraint  are  based  on  the  following  equations. 


nF'F  =TCF'sXi)){RF;s'  +  Rss,s  -  cs  5 rfs) 


S' 5 


iS'SnFS} 


cf’f  ^rc/,'i(o)cs  3c*'(o)r 


iF'  S' 


S'Sr’Sn 


(7.1.11  3) 
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F 


The  time- dependent  terms  in  these  equations  vanish  if  all  of  the  displacements  and  rota 

.F' 

tions  of  S'  relative  to  S  are  along  and  about  6,  (the  axis  of  rotation).  Therefore,  let 


=j?£3f  =  o 

*f f  =<P3F'F  =  0 

The  virtual  displacements  and  rotations  of  the  F'  frame  are  then 


„F 

lF'‘  =  \  0 

l  0 


6RF’F  o  ^  =  ToCf  s\o)  ( 6R$S  +  Ws,SCs'sRfs) 


6Vf - 


F'F 


S’l'p.f 


0 

0 


F'S’n^c.S’S 


T0C*  a(O)«0|, 


(7.1.11-4) 


(7.1.11  5) 
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where  the  use  of  only  the  T0  component  of  T  eliminates  the  time-dependent  terms. 

The  virtual  work  at  5  associated  with  the  force  and  moment  at  F  yields  the  following 
contribution  at  S: 

SW  =(8RSS',S)  Cs'F'(0)T?Fg - 

if  f£+  (7.1.11-6) 

Wi:s)rCs'r'(0)T?M£ 

Dynamic.  The  position  and  orientation  of  the  perturbed  frame  relative  to  the  steady- 
state  position  and  orientation  are  related  as  follows: 


R 


p"  p' 

F" 


-TC 


F"S", 


(0)  (J*S 


F"S 

S" 


+  R 


s"s' 

S" 


(7.1.11-7) 


CF"F'  =rCF"S"(0)  CS"S'  CS'F'(0)  TT 

From  these  equations,  the  virtual  displacement  and  rotation  may  be  obtained.  To  first 
order  in  the  perturbation  quantities, 


si =tcf"s'\o)  -  ~C5'  )<s'] 


(7.1.11-8) 


SrppnF  =TCf  5  (0)  Sxjjgn 


s"s' 


-_s"s' 

where  contributions  due  to  geometric  stiffness  come  from  the  6S„ 

U  is 


term.  The  matrix 


T(f)CF"s  "(0)  -T(t)CF"s"(Q)RF:s' ' 

o  r(t)c,F''s’,(o) 


(7.1.11-9) 


where  the  columns  of  77  correspond  to  J7f,,s  and  0f„s  ,  and  the  rows  correspond  to 


8Rf;f'  and  6^;f' 


Since  77  depends  on  f,  the  time  dependent  terms  must  be  removed  from  the  final 
transformed  equations.  This  is  easily  accomplished  by  taking  the  time- averaged  value 
of  the  transformed  equations.  The  only  contributing  ( i.c .,  nonzero)  terms  then  are  the 
constant  terms,  the  cos2(ftr)  terms,  and  the  sin2(OT)  terms.  In  addition,  since  77  depends 
on  i,  terms  from  matrix  M  will  contribute  to  Jlf,  C,  and  K  in  the  transformed  equations 
and  C  will  contribute  to  C  and  K  by  virtue  of  the  following  relations 


.  p"  p> 

Rpn 


iF"  F' 

0  pn 


+  R 


•.  5"  S' 

RSi, 

•,s"s' 

&S" 


(7.1.11-10) 
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■  ..  F"  F' 

R  p>  i 

■:F"F' 

.  6  pn 


{Qlt 

&s„ 

C 


..  s"  S' 
RSu 

•:s"  s' 

0  C(» 


(7.1.11-11) 


iifG,  the  geometric  stiffness  matrix,  is 


Kg  = 


0  [CS"F"(0  WfFgfRS* 


(7.1.11-12) 


where  the  columns  of  KG  correspond  to  jRj„s  and  0f„s  ,  and  the  rows  correspond  to 
SRf„s  and  6iJ) f"5  . 

7.1.12.  Rotating  Node  Demotion 

The  rotating  node  demotion  constraint  describes  the  relationship  between  two  nodes, 
one  of  which  is  located  in  a  rotating  frame  of  reference  and  the  other  in  a  nonrotating  frame 
of  reference  (fig.  13).  This  constraint  combines  many  of  the  characteristics  of  the  rotating 
frame  and  structural  node  demotion  constraints.  It  is  assumed  that  the  child  frame  F  is 
rotating  about  a  fixed  axis  at  an  angular  speed  fi,  and  that  a  dependent  node  b  is  defined 
relative  to  that  rotating  frame.  The  parent  frame  S  is  stationary  (relative  to  F)  and  an 
independent  node  I  is  defined  relative  to  S.  The  rotating  node  demotion  constraint  is  not 
available  through  the  GRASP  user  interface. 

Steady-State.  The  governing  equations  for  the  steady-state  condition  are  similar  to 
those  for  structural  node  demotion  except  that 


CFS(t)  =T{i)CFS{  0)  =  CF's'(t)  =  T{t)CF's'{  0)  =  CDI{t) 
=r{t)CDI{  0)  =  CD<I\t)  =  T(t)CDI'(0) 

In  addition,  R$s  and  Rf,  1  are  constants,  and 

T(t)  —  To  +  Tc  cos(17t)  +  T,  sin(f It) 


(7.1.12-1) 


(7.1.12-2) 


The  governing  equations  describing  the  deformed  position  and  orientation  of  the  de¬ 
pendent  node  are  then 


-CDF> Cf'fTCfs(0)Css>  [C*I{CII'R$t  +Rjl)  +  R'/  +Rgs}~ 
CDF‘  [Cf'fTCfs(0)Rfs  +  RF‘F  +  R°f'} 

CD'D  =TCD'r{0)CI'ICIS,Cs'sCSF{0)TTCFF'cF,D 


(7.1.12-3) 


In  order  to  be  independent  of  t,  let  T  —  To  and  choose  Rp,F'  =  0  and  CDF‘  —  A.  Thus, 
only  displacements  along  and  rotations  of  D  about  the  axis  of  rotation  can  be  nonzero. 
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Figure  13.  Rotating  node  demotion  constraint. 


The  virtual  work  at  D'  due  to  a  force  and  a  moment  acting  at  D'  is  determined  in  terms 
of  the  virtual  displacement  and  rotation 


«ig'u  =  -  CDF'Sfl',FC’”r%CFS(0)Css'  +  CllF'CF'FTCFS{0)Css'R%.s'6j,s,S± 
cnr" CF' FT<,CF'\0)CS“’  [cs' '(ty'j ' c"' Rfi' ''  +  SR',1')  + 

CDF'{  -  &irp,FCF'FT,CFS(0)R^s  + 


=cPFT,crs(o)cs,ti’j '  +  c‘"'r0c''s(o)c5s'6^|;5  -  curtyfe 


i  FSr 


UFm  nFS< 


iSS'  cj.  S' S 


*DF'i„lF'F 


(7.1.12-4) 
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The  virtual  work  done  at  D'  is  then 


6W  -  {5RFF',Ff  Fj?:  -  {S1>£FY  {Mpl  +  R°:fF°')+ 
(6Rss',sjrCs'sCSF{0)T0Fj?'+ 

(67l)f,s)r[Cs  sCSF(0)TF Mp'  +  kg', s Cs' s Csf(0)Tq  Fp' ]+  (7.1.12-5) 

(^f;)Tc/scSF(o)r0F^'+ 

(6ti'I)T  [ClsCFS(Q)T0rM}?'  +  Rj*'1' CISCSF(0)TFFp) 

Dynamic.  The  governing  equations  are  similar  to  those  for  structural  node  demotion 
except  that,  as  in  the  static  case,  RF,F  =  0  and  CDF  =  A.  The  governing  equation  for 
the  position  of  the  perturbed  dependent  node  is 


where 


rd"d'  =RD"r  +  Err  +  Rn  +  Rrs"  +  Rs-s'  +  Rs>s_ 

rns  _  rfn  _  RF'F  _  rf"  F  _  rdf"  _  RD'D 


rns  +  rfn  ^  RS'S  +  RN'S>  +  RF'N'  _  RF'F 


Solving  for  the  dependent  node  displacement, 

RD"D'  =rU"1"  +  +  R1'1  +  RlS"  +  RS"S '  _  rN'S'  _ 

RF'”'  -rF"f'  _rDF"  _rd'd 

By  referring  the  displacements  to  the  D  basis  this  vector  relation  becomes 

rV"d‘  =CDF"cF''F'cF'N'CN's'Cs,s"cs"ICIICI'I"RfnI"  + 

CDF"CF"F‘CF'N'CN'S'CS'S"CS"I{RI"I>  f  RI'I')+ 

cdf"cf"f-cf>n'cn-s'cs>s"{riss"  Rs"s')- 

gDF"  £F"F'  qF'N' qN'S'  Rn'S'  _ 

CDF"  CF"F'  CF'N'  RFjN'  _  CDF"  (Rp"F'  +  R<j>F"  )  -  R%D 
The  virtual  displacement  is  then 

6Rg'D'  =CDF,'CF',F,CF'N,CN,S'CS'S''CS''I6RI;,I+ 

cDF"cF"F-cF-N'cN'S'cS'S"6RSy>  _  CDF" SRFp„F' + 
C^F" CF" F1  CF' N‘  CN' S'  CS‘ S"  CS" I CII'  1" RD"l" 

CDF''CF"F'CF'N'CN'S'6CS"S'{CS"ICJI'CI'I "Rffi"  + 

CDF"6CF''F’CF'N'{CN'S'CS'S''CS''ICII'C1'I "rd;;i"  + 
cJV5'cS*S"c5»/(ji/"1'+jR/*/)  + 

C^S'CS'S"{RISS"  +  rS"s')  _  cn,5'r^:s'  -  Rf'”'} 


(7.1.12-6) 


(7.1.12-7) 


(7.1.12-8) 


(7.1.12-9) 


(7.1.12-10) 
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and  the  virtual  work  done  by  the  force  at  D'  is 


(, srz"d')  f% '  hsri"1')  iciN'cN,F'Fj?; + ciN\cN,F'Fp:ycN's’ef:s'- 

rilN'nN'F'  riD'pF"F'\  , 


CIN  CN  F  Fg,  6£„F  J+ 

(6rss:s')t{cs'n,cn'f,fp: + CS' N\CN'F'F?:y  CN'S' elv¬ 
es' N'  cN'F'  fD'qF'i'F']_ 

{srf;:f,)tf^:-v 

ciN,RD:F(cN'F‘Fj?:)-cN's'oss:s'  -  cin'rd7cn'f,f?;9f:f')+ 
(Sipfn  s>  )T[Cs  n> Rtf!1'  CN  F'  Fp!  -  Cs'n'(Cn'f'fP:)'Cn'i'r1j"1'~ 

cs'n,(cn,f'f°;)~cn's,rss:s,+ 

cs'N,(cN'F'Fp;ycN,iR?,i'op,+ 

cs'N'R%i'(cN'F,Fp:ycN's,dss:s'- 

(8i>£:F')Ti(cF,”'R°,F'yFp;  -  p^c^c^r^- 

f°:cf,n'cn,s'rss:s' + Fp;cF'N'cN,iR,j>,i,9i]"i,+ 

pD'  nF'N'  rN'S‘  f>D'S'pS"S'  £ D'(nF'N 1  nD' F1  pF"  F'  1 
r  pi  O  O  llgi  Ugii  —  JF  pi  It.  j\ji  )  Opn  | 

(7.1.12-11) 

The  virtual  rotations  can  be  easily  obtained  from  the  variation  of  the  direction  cosine 
relation 


8Cd"d'  =  -  CD"D,CD'D8iD  CDD' 

=CD''/,'5C/',/'c/'/C/s''cs''5'cs'N'cN'F'c,:’'F''ci:’''DCDD,-|-  (7.1.12-12) 

qD" I"  gl" /'  Ql' I gis" £(75” S'  (jS' N'  f-iN' F'  qF' F"  qF"  D gDD'  + 
qD"1"  qI'IqIS"  (jS"S‘  qS'N'  qN'F'  fiCF'F"  qF"  DqDD' 


thus  yielding 


_ i"  r' 

iDlcj,1  1  nlD  ,  nDS 


64'D  =  Cu16tP)  ’  CIU  +  CVS'MW  cS"D  -  CDF'l6ippn '  c 

Upon  removal  of  the  tilde  the  virtual  rotations  are 


Wg"D'  =CDF" CF"F' CF'N'C”  *  C*  5” C*  ' Hi  1  + 

cDF''cF''F'cF'N'cN'S'CS'S''HSy  _  CDF"  b^'f 


(7.1.12-13) 


(7.1.12-14) 
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and  the  virtual  work  done  by  the  moment  is 


wg'D')TMg  a  -  9|;;s')cs'w'cw'f'(a + 0£:f')m?;+ 

(Si,§;:s'f(  a  -  0§::s‘)cs'n'cn'f(a + f;:,f')m%  - 

=(Si,‘,'‘'')T\cw'ctl'F'M$  +  c>N\cN‘F'M?:)'cK's'ef:s'- 

C'n,Cn'f'm$6f:f'\+ 

W§:s')t{cs'n'cn'f'm!?;  +  cs'N’(cN'F'M^)'cK's'e 
cs’N'cN'F'M^;ef,F'i- 
(»£:F')TMg 

(7.1.12-15) 

since  CF  N  =  T,  let  CN  5  =  CF  s  ,  as  in  the  structural  node  demotion  constraint.  The 
time-dependent  terms  can  be  eliminated  when  RpF  =  0  and  CDF  =  A. 

Combining  the  virtual  work  due  to  the  force  and  the  moment  at  P",  the  matrix  % 
can  then  be  calculated. 


n 


qDF'j<qF'i  -Cdf'TCf'iRY'1'  Cdf>TCf's' 
0  Cdf'TCf>i  0 

—Cdf'tCf's' Rg'i s'  -CDF'  0 

CDF'TCF'S'  0  _CDF' 


(7.1.12-16) 


where  the  columns  of  71  are  associated  with  SRj"1  ,  St/)}"1' ,  SR§”s\  6ifi§',',s' ,  8Rp',F,  and 
8ipF„F  ,  respectively;  and  the  rows  correspond  to  SRp  D  and  Sip®  D  ,  The  coefficients  of 
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the  matrix  KG  are 


6R 


rr 

I 


T 


row: 


row: 


R\  1  column: 
6j  1  column: 
RgifS  column: 

djl#  nl 

dg„  column: 

pH  p> 

Rpt,  column: 
6p„ F  column: 

ill  ji 

Rj  column: 
61,"1'  column: 
R%,S  column: 
0g„ s  column: 
RpnF  column: 


0 

0 


u  (7.1.12-17a) 

0 

cif'ttfP,' 


0 


c/F'(rTF^')'cF,/ 


(7.1.12-17b) 


[ rYvC,f\TtF$)  +  C,I‘"{T1Mp 


,IF‘t 


D' 


)  )C* 


F'  S' 


0 


0f:f'  column:  R? V ClF> Tt F$  +  Cif’TtM$ 
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qII  c1  7 

SRss„s 


row: 


Ri  1  column:  0 


&\  1  column:  0 


Rss„s  column:  0 


ess:s'  column:  -  CS'F’ (Tt F$)  CF’ s' 


Rp„F  column:  0 


eF;:r  column:  Cs'f'TtFP: 


F' 


Sip 


s"  s' 

S" 


row: 


R\  r  column:  CSF(TTFg) C 


s' F'  (rpT  r*D' \~  r'F' I 


e1,"1'  column:  -  CS'F' {T* F$)  CF’ T 


Rss:s'  column:  Cs' F' (Tt  Fg )' CF'S' 


0ss:s'  column:  -  [R$s‘ CS'F' (TT Fg )  +  CS  /'{T 


nM  17*.  . 

Rp„  column:  0 


6f:f'  column:  R$s' CS'F' Tt  Fg  +  Cs' F' Tt  tig 


(7.1. 12-1 7c) 


f7.1.12-17d) 
Alg  )~}CFS' 
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row: 


i2|,/S  column:  0 
<9 J„  5  column:  0 
Rp„F  column:  0 


(7.1.12-17e) 


(7.1.12— 17f) 


0F„F  column:  0 

7.1.13.  Rotating  Generalised  Coordinates 

The  rotating  generalized  coordinate  constraint  relates  generalized  "oordinates  in  one 
subsystem  to  the  corresponding  generalized  coordinates  in  another  subsystem  that  is  ro 
tating  at  constant  angular  speed  relative  to  the  first.  This  constraint  is  often  applied  to 
subsystems  that  contain  periodic  structures.  This  constraint  is  not  available  through  the 
GRASP  user  interface. 

Steady- State.  The  general  form  of  the  transformation  from  rotating  to  nonrotating 
coordinates  is 

<M  =  Tqw  (7.1.13-1) 

where 

T  =  T0  l-  Tc  cos ft*  -|-  T,  shift*  (7.1.13-2) 
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and  qn  is  a  set  of  generalized  coordinates  in  the  rotating  subsystem  that  corresponds  to  a 
set  of  qn  generalized  coordinates  in  the  nonrotating  subsystem. 

In  order  to  make  the  transformation  equations  independent  of  time,  let  T  =  To  and 
SqRT  —  TySq jv>  This  eliminates  any  generalized  forces  of  the  lateral  (cosine  or  sine)  type. 
Then,  the  virtual  work  is 

&W  =  qlQn  =  qJrT?QR  (7.1.13-3) 

Dynamic.  In  the  rotating  system,  the  virtual  work  can  be  written  as 

5W  =  qR[MRqR  +  CRqR  +  KRqR)  (7.1.13-4) 

where 

qn  =TqN 

qR  -TqN  4-  TqN  (7.1.13-5) 

qR  -TqN  +  2  fqN  +  TqN 

Substituting  these  relations  into  equation  (7.1.13-4),  the  virtual  work  can  be  obtained  in 
terms  of  the  generalized  coordinates  of  the  nonrotating  system. 

—  ‘In'T'1'  MRTqft  +  (CrT  -f  2MRT)qN  4-  (KRT  4-  CRT  -f  MRT)qw  (7.1.13-6) 

Thus,  the  C/v  coefficient  matrix  (in  the  nonrotating  system)  depends  on  CR  and  MR,  and 
Kn  depends  on  Kr ,  CR,  and  Mr. 

7.1.14.  Rotating  Air  Mass 

The  rotating  air  mass  constraint  transforms  the  air  node  generalized  coordinates  and 
their  associated  generalized  forces  between  a  rotating  subsystem  and  a  nonrotating  sub 
system.  As  in  the  other  air  mass  constraints,  the  air  node  generalized  coordinates  are  not 
transformed  out  of  the  inertial  frame  of  reference.  The  rotating  air  mass  constraint  is  not 
available  throught  the  user  interface  in  GRASP. 

Steady-State,  oincc  only  Uf  and  ifr  are  valid  coordinates  in  the  steady-state  prob¬ 
lem,  and  both  are  rotationally  symmetric,  they  are  treated  in  exactly  the  same  manner  as 
in  the  copy  air  mass  constraint. 

Dynamic.  For  a  set  of  dynamically  perturbed  air  node  generalized  coordinates,  let 

f  f A  1  f  «A  ) 

qn  =  S  ?  i  S<1R  =  )  W12 r  ?  (7.1.14-1) 

l.  0J3R  j  1  ^13 R  j 
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7.2.  Composite  Constraints 

In  general,  a  composite  constraint  is  a  constraint  that  is  built  up  out  of  one  or  more  of 
the  primitive  constraints  that  have  been  described  in  the  previous  sections.  The  bundling 
of  primitive  constraints  into  a  single  constraint  is  primarily  done  for  the  convenience  of  the 
user.  There  are.  many  times  that  sets  of  constraints  must  be  used  together,  and  it  mokes 
sense  to  combine  them  internally.  In  the  following  sections,  the  composite  constraints  that 
have  been  constructed  from  the  set  of  primitive  constraints  in  GRASP  will  be  described. 
All  of  the  composite  constraints  are  available  from  the  GRASP  user  interface. 

7.2.1.  Aeroelastic  Beam  Connectivity 

The  purpose  of  the  aeroelastic  beam  connectivity  constraint  is  to  provide  a  means  for 
attaching  an  aeroelastic  beam  element  to  a  GRASP  model.  The  element  subsystem  for 
the  aeroelastic  beam  consists  of  a  frame  of  reference,  a  root  node,  a  tip  node,  and  an  air 
node,  all  of  which  must  be  connected  to  their  counterparts  in  the  existing  portion  of  the 
model.  Therefore,  the  aeroelastic  beam  connectivity  constraint  must  contain  a  fixed  frame 
constraint  (for  the  frame),  two  structural  node  demotion  constraints  (for  the  root  and  tip 
nodes),  and  a  copy  air  mass  constraint  (for  the  air  node). 

In  the  definition  of  the  aeroelastic  beam  connectivity  constraint,  the  position  ami 
orientation  of  the  dependent,  element  root  node  i?  relative  to  an  existing,  independent  node 
In  must  be  provided.  The  position  and  orientation  of  In  relative  to  its  subsystem  frame  of 
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reference  (the  superordinate  frame  Sr)  is  known  from  the  definition  of  Ir.  Therefore,  the 
position  and  orientation  of  the  dependent,  root  node  relative  to  the  superordinate  frame 
can  be  calculated. 

qRSr  „qRIr(jIrSr 


(7.2. 1-1) 

R™*  -CSnInR^R  + 


After  locating  the  parent  subsystem  of  the  element  subsystem  in  the  system  organization 
tree,  the  position  and  orientation  of  the  parent  frame  relative  to  the  superordinate  frame 
can  be  calculated.  Since  the  element  frame  and  the  element  root  node  are  coincident, 
the  position  and  orientation  of  the  element  frame  relative  to  its  parent  frame  can  then  be 
determined. 


qPP  qRSr  qSrP 


(7.2.1-2) 

Rfpp  =RP'  =  +  Rssf) 


With  this  information,  the  fixed- frame  constraint  can  be  defined.  In  addition,  all  of  the 
position  and  orientation  information  is  available  to  define  the  structural  node  demotion 
constraint  for  the  clement  root  node.  In  those  cases  where  the  superordinate  frame  is 
not  the  same  as  the  parent  frame,  it  is  necessary  to  create  copies  of  the  independent  and 
element  root  nodes  in  each  of  the  subsystems  leading  to  their  nearest  common  ancestor. 
These  nodes  are  chained  together  by  a  series  of  structural  node  demotion  constraints. 

The  position  of  the  clement  tip  node  T  relative  to  the  root  node  is  defined  as  R1  n  - 

ib3  ,  and  the  orientation  Crn  is  defined  as  an  Euler  rotntion  of  magnitude  8'(  about  . 
After  the  position  and  orientation  of  the  root  node  relative  to  the  independent  tip  node 
1 1  has  been  calculated,  the  offset  of  the  element  tip  nude  from  the  independent  tip  node 
can  be  determined. 


RTh  ==c/rSr(CSTni?T«  +  RRST  __  r!tSt ) 
qTIt  qTRqRSt  qSt  It 


(7.2.1-3) 


At  this  point,  the  structural  node  demotion  constraint  for  the  element  tip  node  can  be 
defined.  In  those  cases  where  the  superordinate  frame  is  not  the  same  ns  the  parent 
frame,  it  is  necess«.»y  to  create  copies  of  the  independent  and  clement  tip  nodes  in  each  of 
the  subsystems  leading  to  their  nearest  common  ancestor.  This  creates  another  chain  of 
structural  nodes,  all  connected  together  by  structural  node  demotion  constraints. 

If  the  beam  element  is  to  be  connected  to  an  uir  mass  element,  the  position  and 
orientation  of  the  clement  subsystem  relative  to  the  corresponding  air  nude  is  calculated. 
Then,  the  copy  air  mass  constraint  is  defined.  In  those  cases  where  the  air  node  is  nut 
defined  in  the  patent  frame,  it  is  necessary  to  create  copies  of  the  independent  and  element 
air  nodes  in  each  of  the  subsystems  leading  to  their  nearest  common  ancestor.  This  creates 
a  chain  of  air  nodes,  all  connected  together  by  copy  air  inass  constraints.  If,  however,  the 
beam  clement  is  not  to  be  connected  to  an  air  mass  element,  the  four  nodal  air  mass 
degrees  of  freedom  are  constrained  out  using  prescribed  constraints. 
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7.2.2.  Air  Mass  Connectivity 

The  purpose  of  the  air  mass  connectivity  constraint  is  to  provide  a  means  for  attaching 
an  air  mass  element  tu  a  GRASP  model.  The  air  mass  element  subsystem  is  unusual  in 
that  the  frame  serves  only  to  establish  the  position  and  orientation  of  the  element  relative 
to  the  remainder  of  the  model.  Therefore,  while  the  frame  does  exist  and  docs  need  a 
frame  constraint,  it  has  no  frame  degrees  of  freedom.  The  air  mass  connectivity  constraint 
is  then  made  up  of  a  fixed-frame  constraint,  a  copy  air  mass  constraint,  and  one  or  more 
prescribed  constraints. 

In  the  air  mass  connectivity  constraint  it  is  assumed  that  the  independent  air  node  /, 
the  dependent  (clement)  air  node  A,  and  the  element  frame  F  are  all  coincident. 


CIA  -CAF  -  A 
R,a  =Raf  =  0 


(7.2.2-1) 


After  locating  the  parent  of  the  element  subsystem,  the  position  and  orientation  of  the 
pnicnt  frame  relative  t.o  the  superurdinate  frame  can  be  calculated.  Since  the  position  and 
orientation  of  the  independent  air  node  relative  to  its  subsystem  frame  (the  superordinate 
frame  5)  is  known,  the  position  and  orientation  of  the  element  frame  relative  to  the  parent 
frame  can  also  be  calculated. 


CFF  -cip  -  C,SCSP 
R?.p  =R‘PP  »  CPS(R1sS  +  Rf) 


(7. 2. 2-2) 


These  expressions  provide  the  information  necessary  to  define  the  fixed  frame  con 
strnint.  In  addition,  the  copy  air  mass  constraint  can  be  defined  at  this  time.  In  those 
cases  where  the  superurdinate  frame  is  not  the  same  ns  the  parent  frame,  it  is  necessary  to 
create  copies  of  the  independent  and  element  air  nodes  in  each  of  the  subsystems  leading 
to  their  nearest  common  ancestor.  These  additional  air  nodes  are  also  connected  together 
using  copy  air  mass  constraints. 

If  the  model  containing  the  air  mass  connectivity  constraint  is  to  be  used  in  n  steady 
state  problem,  the  two  cyclic  aii  nude  degrees  of  freedom  aie  meaningless.  Therefore,  they 
must  be  eliminated  by  defining  two  prescribed  constraints  in  the  supcrordinati  subsystem. 
If,  on  the  other  hand,  the  model  is  to  be  used  in  an  eigensolution,  the  gradient  degree  of 
freedom  is  meaningless.  A  single  prescribed  constraint  is  then  defined  in  the  supc lotdinaU 
subsystem. 
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7.2.3.  Periodic  Structure 


The  purpose  of  the  periodic  structure  constraint  is  to  provide  a  simple  means  for 
creating  an  axially  symmetric  structure.  This  is  accomplished  by  replicating  a  single  branch 
of  the  model  at  equal  azimuth  angles  about  an  axis  of  symmetry.  For  this  constraint,  the 
parent  subsystem  represents  the  assembled  periodic  structure  and  the  child  subsystem 
represents  a  single  component.  The  periodic  structure  constraint  consists  of  one  or  more 
of  the  following:  a  periodic  node  demotion  constraint,  a  periodic  generalized  coordinate 
constraint,  and  a  periodic  air  mass  constraint.  Note  that  the  periodic  frame  constraint 
must  be  defined  separately. 

When  there  are  nodes  in  the  component,  periodic  node  demotion  constraints  are 
needed  to  transform  them  into  the  assembled  structure.  If  the  independent  node  corre¬ 
sponding  to  a  dependent  node  (in  the  component)  does  not  exist  in  the  parent  subsystem, 
a  string  of  images  of  the  independent  node  are  created  in  the  intervening  subsystems  and 
chained  together  with  structural  node  demotion  constraints.  Similarly,  if  the  dependent 
node  does  not  exist  in  the  child  subsystem,  a  string  of  images  of  that  node  are  created 
and  chained  together.  Since  the  independent  node  (or  its  image)  now  exists  in  the  parent 
subsystem  and  the  dependent  node  (or  its  image)  exists  in  the  child  subsystem,  a  periodic 
node  demotion  constraint  can  be  defined. 

One  or  more  periodic  generalized  coordinate  constraints  are  needed  if  there  are  gen¬ 
eralized  coordinates  in  the  child  subsystem.  Similarly,  one  or  more  periodic  air  mass 
constraints  arc  needed  if  there  arc  air  nodes  in  the  component.  A  process  identical  to  that 
used  to  connect  structural  nodes  is  used  if  the  dependent  and  independent  air  nodes  are 
not  in  the  child  and  parent  subsystems,  respectively. 

7.2.4.  Rigid-body  Connection 

The  purpose  of  the  rigid-body  connection  constraint  is  to  provide  a  simple  means  for 
connecting  two  nodes  together  rigidly.  It  is  actually  a  special  case  of  the  screw  constraint 
in  which  the  translation  and  rotation  degrees  of  freedom  are  both  locked. 

7.2.5.  Rigid-body  Mass  Connectivity 

The  purpose  of  the  rigid-body  mass  connectivity  constraint  is  to  provide  a  means 
for  attaching  a  rigid-body  mass  element  to  a  GRASP  model.  The  clement  subsystem 
consists  of  a  frame  oi  reference  and  a  center-of-mass  node,  both  of  which  must  be  connected 
to  their  counterparts  in  the  existing  portion  of  the  model.  Therefore,  the  rigid  body 
mass  connectivity  constraint  is  made  up  of  a  fixed  frame  constraint  and  a  structural  node 
demotion  constraint. 

In  the  definition  of  the  rigid-body  mass  connectivity  constraint,  the  position  and 
orientation  of  the  dependent,  element  center-of-mass  node  C  relative  to  an  independent, 
existing  node  I  is  provided.  The  position  and  orientation  of  the  independent  node  relative 
to  its  subsystem  frame  of  reference  (the  superordinate  frame  5)  is  also  known.  Then,  the 
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position  and  orientation  of  the  center- of- mass  node  relative  to  the  superordinate  frame  can 
be  written. 


£jCS _ qCI  qI  s 


Rf s  =CSIRf1  +  4s 


(7.2. 5-1) 


After  locating  the  parent  subsystem  of  the  element  subsystem  in  the  system  organization 
tree,  the  position  and  orientation  of  the  parent  frame  relative  to  the  superordinate  frame 
can  be  calculated.  Since  the  element  frame  and  the  element  center-of-mass  node  are 
coincident,  the  position  and  orientation  of  the  element  frame  relative  to  its  parent  frame 
can  then  be  determined. 


CFP  _ qCP  _  qCSqSP 

Rpp  =RCpP  =  Cps{Rcss  -1-  Rssp ) 


(7. 2. 5-2) 


With  this  information,  the  fixed-frame  constraint  can  be  defined.  In  addition,  all  of  the 
position  and  orientation  information  is  available  to  define  the  structural  node  demotion 
constraint  for  the  center-of-mass  node.  In  those  cases  where  the  superordinate  frame  is 
not  the  same  as  the  parent  frame,  it  is  necessary  to  create  copies  of  the  independent  and 
element  ccntcr-of-mass  nodes  in  each  of  the  subsystems  leading  to  their  nearest  common 
ancestor.  These  nodes  are  chained  together  using  a  series  of  structural  node  demotion 
constraints. 


7.2.6.  Rotating  Structure 

The  purpose  of  the  rotating  structure  constraint  is  to  provide  a  simple  means  for  allow¬ 
ing  one  subsystem  to  rotate  relative  to  another.  For  this  constraint,  the  parent  subsystem 
represents  the  nonrotating  structure,  while  the  child  subsystem  represents  a  rotating  struc¬ 
ture.  The  rotating  structure  constraint  consists  of  one  or  more  of  the  following:  a  rotating 
node  demotion  constraint,  a  rotating  generalized  coordinate  constraint,  and  a  rotating  air 
mass  constraint.  Note  that  the  rotating  frame  constraint  must  be  defined  separately. 

When  there  are  nodes  in  the  rotating  subsystem,  rotating  node  demotion  constraints 
are  needed  to  transform  them  into  the  nonrotating  subsystem.  If  the  independent  node 
corresponding  to  a  dependent  node  (in  the  rotating  subsystem)  does  not  exist  in  the  parent 
subsystem,  a  string  of  images  of  the  independent  node  are  created  in  the  intervening 
subsystems  and  chained  together  with  structural  node  demotion  constraints.  Similarly,  if 
the  dependent  node  does  not  exist  in  the  rotating  subsystem,  a  string  of  images  of  that 
node  are  created  and  chained  together.  Since  the  independent  node  (or  its  image)  now 
exists  in  the  parent  subsystem  and  the  dependent  node  (or  its  image)  exists  in  the  child 
subsystem,  a  rotating  node  demotion  constraint  can  be  defined. 

One  or  more  rotating  generalized  coordinate  constraints  are  needed  if  there  are  gen 
eralized  coordinates  in  the  child  subsystem.  Similarly,  one  or  more  rotating  air  mass 
constraints  .ue  needed  if  there  aic  air  nodes  in  the  component.  A  process  identical  to  that 
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used  to  connect  structural  nodes  is  used  if  the  dependent  and  independent  air  nodes  are 
not  in  the  child  and  parent  subsystems,  respectively. 
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8.  ELEMENTS 


The  GRASP  element  library  currently  contains  three  elements,  the  aeroelastic  beam, 
the  air  mass,  and  the  rigid-body  mass. 


8.1.  Rigid-Body  Mass 

In  GRASP,  rigid  bodies  are  modeled  as  being  influenced  only  by  inertial  and  gravita¬ 
tional  forces. 

For  the  purposes  of  modeling  the  motion  of  a  rigid  body  in  an  inertial  and  (possibly) 
gravitational  field,  consider  a  rigid-body  mass  element  D  tlmt  has  an  inertia  dyadic  L 

^teady-State.  The  rigid-body  mass  element  (fig.  14)  has  a  body-fixed  node  N  and  a 
frame  of  reference  F.  Node  N  is  initially  coincident  with  the  deformed  frame  F'  ( RF  N  —  0 
and  CF  N  —  A).  The  virtual  work  at  the  deformed  node  N1  is 

6W  =  Fn‘  •  iSRn‘i  +  MN<  •  fy"’1  (8.1-1) 


from  which  nodal  forces  and  moments  can  be  derived.  The  nodal  virtual  displacment  and 
rotation  variables  for  this  element  are  SRfl  N  and  6ip$  N ,  respectively. 


The  inertial  virtual  displacement  and  rotation  of  the  deformed  node  N1  are 
lbRN>I  =  W(/  +  HF>I  X  Rn  n  +  n6Rn'n 


5i,N,I=S^'N+6^NF'+6rPF'1 


(8.1-2) 


The  force  acting  on  the  body  at  N1  is 

Fn>  =  -mAN>I  +mg  (  (8.1-3) 

where  the  inertial  acceleration  of  N'  is 

an'‘  =  af''  +  nF''  x  (nF'1  x  (8.1-4) 

Substituting  equation  (8.1  4)  into  equation  (8.1  3),  and  transforming  from  the  body  fixed 
(Ar)  coordinate  system  into  the  dcformed-frame  (F1)  system, 


Fpf  =m{gN  -  A^1  -  R%N) 

iF'I  nF'InF'InN'N') 


Z  I?1  T  ~  V*1 


(8.1-5) 


■  m 


br.-AP-nFn£lRXN) 
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Figure  14.  Rigid-body  mass  element. 


The  moment  acting  on  the  body  at  N1  is 

mn'  =  hn>i  x  aN>I 

where  the  inertial  angular  momentum  at  N'  is 

h"'1  =  I- n"'1 

and  the  inertial  angular  velocity  at  N'  is 

n"'1  =  nN'N  +  nNF'  -m?'1 


Substituting, 

where  lift1  =  . 


Mg' 


-  ijN'i 

—  it  pi 


n 


N'l 

F< 


The  frame  force  and  moment  components  can  be  derived  in  a  similar  manner.  If  the 
frame  virtual  displacements  and  rotations  are  6RFF  and  SipFiF ,  respectively, 


jpF'  _  r<F'N  jpN' 
Jc  pi  —  O  jTm 


N 


-  7?W 


Mp!  =  Mg'  +  R%'NF$' 


frN'  tjN'N 
-  tN  liN 


(8.1-10) 


Dynamics.  Assuming  that  the  rigid-body  mass  node  is  perturbed  from  its  steady- 
state  position  (fig.  14),  the  virtual  work  at  the  node  may  be  expressed  as 


SW  =  6R%’>TF»"  + 

(8.1-11) 

where  the  force  and  moment  are 

F$'  =  -  mA^"1  +  mgN 

(8.1-12) 

Mf  =  -  IjvaJJ"'  +  H8UtSlg'1 

the  angular  momentum  is 

h!T  =  /*<' 

(8.1-13) 

and  the  inertial  angular  velocity  is 

nN'7  _  o F"F'  i  nN"N'  ,  rF"F'c>F'I 

~r  ifp’i 

(8.1-14) 

The  virtual  displacement  is  then 

ajijp  =  srff:f'  +  +  Hp^' [rff:f'  +  r%"n  +  rnn"} 

and  the  virtual  rotation  is 

S^"1  =  6^'  +6^" N' 


(8.1-15) 


(8.1-16) 
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The  inertial  acceleration  of  the  node  is 


7  =Ap'1  +n£tsi£tR%1'  +  rf„ 


..  F"  F‘ 


~  -i N"N‘ 

-  R&  ’'Or,,  +  R, 


WN(1FF:'  +  + 


\F'InN'N 


N  VF" 

t  l?tl  pi 

~ F'lhN'N -yX 


lN 


4"  2fljr  i  7  Rpn  — 


,  F"F' 


^f'i  f}N"N'  +  rf;f' + 


2 Rpj 

iaff:j  -  (npRg")'^1  - 

ng'ng'Ajr1"* 

and  the  components  of  the  gravitational  acceleration  are 

^  »  nil  nl 

9N  =  9F'  +  9F'9F" 


(8.1-17) 

(8.1-18) 


The  components  of  the  inertia  dyadic  in  the  nodal  basis  can  be  expressed  in  matrix 
form  as 


-N"N‘  -N"N' 

In  =  (A  +  eN  )CNN'lN„CN'N{  A-6n  ) 


(8.1-19) 


Finally,  the  angular  acceleration  is 


„  ,  ■■F"F'  zN^N' 


a 


N  N 
N 


-i .  iF"F'  _  ,  iN"N' 
F  In  i  nF  Ir 


9  pn  4*  4*  &  pi  0  pu  +  nF,  9 


N 


(8.1-20) 


The  force  and  moment  can  then  be  obtained  from  the  substitution  of  equations  (8.1- 
12)  through  (8.1-20)  into  (8.1-11).  When  the  virtual  work  is  calculated,  it  consists  of  the 
same  steady-state  residuals  {(?}  as  were  obtained  in  the  previous  section,  in  addition  to 
the  virtual  work  associated  with  the  coefficient  matrices  [M],  [C],  and  [if]. 


(srf:f' 

I 

{  6R""N' 

Wr*' 


'  ••  F"  F'  ' 

'  •  F"  F'  ' 

Rpn 

-F"F' 

R  Fn 
.  F"  F' 

(  nF"  F'  N 

ltFn 

l*f) 1 

9f„ 

•:  N"N' 

Rn 

•iN"N' 

1  +  [0]  < 

0F„ 

■.N"N' 

rn 

■„N"N' 

>  +  [if]  < 

d£F' 

R%“N' 

nN"N' 

UN  ; 

Jn 

4 

(8.1-21) 


78 


where  the  M  coefficient  matrix  is  defined  to  be 


SRl 


F" 


row: 


..  f"p' 

Rpn  column:  m 


..p"  p1 

9Fi,  column:  —  mR 


-N'N 

N 


(8.1-22a) 


c  /  F  F 
vippn  row: 


•i  N"N' 

Rn  column:  m 


"N"  N 

0N  column:  0 


•i  F"F‘  .  , 

Rp„  column:  miZjv 


CF  column:  -  mR^'NR^N  +  CNN'lN..CN'N 


(8.1-22b) 


•:  N"N'  , 

Rn  column:  mR^ 


■iN"n 

dN  column:  CNN  IN»CN  N 


6R%"N'  row: 


..  p"  F' 

Rp„  column:  m 


»F"F' 


0  pa  column:  - mil $ 


N'N 


(8.1-22c) 


■iN" N ' 

Rn  column:  m 


•:N"N 

0N  column:  0 
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column:  0 


Sift 


N 


row: 


•i  F"F' 
Rpu 


fj  column:  CNN'lN»CN'N 


«  N"N' 

Rn  column:  0 


9n  column:  CNN  In<iCn  N 
tlic  C  coefficient  matrix  is  defined  to  be 


SRF„F  row: 


F"  F' 


F'l 


Rpu  column:  2m$lp, 


*  F* 

0F„  column:  —  2mflp,IR^N 


■„N"N 


F'l 


Rn  column:  2  mwF, 


iN"N 

6n  column:  0 


Siftp„F  row: 


.p"p 


Rpu  column: 


.  p"  p' 


'eF„‘  column:  -  2mR%'"np R%’N  -  (CNN' IN..CN'Ntt %' 

cNN'iN..cN’Nn$,‘  +  a%lcNfl'iN.,cN‘N 


N"N' 

Rn  column: 


lN"N 

0N  column: 


2mR^i'N(lp,1 

U1plICNN'lN,<CN'N  ~ 


(8.1-22d) 


(8.1-23a) 


)  + 

(8.1-23b) 
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row: 


SR»"n' 


•„F"F' 

Rpu 

column: 

2  mflp,1 

,F"  p' 

-2  mhFF:rR^N 

0  pn 

coluiim: 

•_N"N: 

Rn 

1 

column: 

2mflp/ 

•_NnN 

8i>%  N  row: 


8n  column:  0 


.  p"  p' 

Rp„  column:  0 


I F  column:  CN!f'lN,,CN'Knp +{l^CNN'lN.,CK'N- 

(cKN'iN„cN'Nnf,,y 


lN"N' 

Rn  column:  0 


■  N"N 


eN  column:  in..cn'n  -  (cNN' iNl.cN'Nn%ry 

and  the  K  coefficient  matrix  is  defined  to  be 


r  r%Pn  F* 

oRpn 1  row: 


-  -  ...  -  *v  nf  r  •  nl  r 

Rp„  column:  mUF,  Up, 

9p„l‘  column:  mlAp,1  —  gp>  -  (Up,1  Rp/ N)  Up,1- 

ttpR%'Nn%t) 

N  column:  mUp^Up,1 
tip,  K  column:  0 


(8.1-23c) 


(8.1-23d) 


+ 


(8.1-24a) 
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row: 


Sip 


pn 


Rp,F  column:  —  mlAp,1  -  §f'  +  (ftjr/fip/.Rjv  A)  - 

AJ5,Jvnf:jrnf:xi 


dF„F  column:  mR^  N[AFiI  -  gF1  —  N)  ftp'1— 

a&'sg"  ft*;.1)  +  nZ'cw'ijv.'C^np- 

(CWN,IJV«Cw'Nfl?)'flp 


8R 


N 


Rfti  N  column:  —  —  gF1  4-  (ftp^ftp^R^f  N) 

rX"^1^1) 

ofN'  column:  n&CNN‘ls»CN'Nil£I- 

n£'(c™Vcw,iVnST 


row: 


Rf„f  column:  mftp/ftp. 


-  F'lnF'I 


0F"F'  column:  rr^Ap,1  -  §F'  ~  {ftp^R^  N)  ftp.1- 

n(,IRglta&1) 


Rh  N  column:  mftp/ftp,1 


6tJ>m  n  row: 


0$  N  column:  0 


Rp"F'  column:  0 


6f:F‘  column:  ft^C"”' - 


t pi 

iNN1 


(C™  IN..C"  ”ft£'1)  ft 


F< 


Rp,  N  column:  0 


0»"N  column:  ft^C^'l^^ft^1- 


6  pi 


F> 

N’NnF'I) 


(8.1-24b) 


(8.1-24c) 


(8.1-24(1) 
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8.2.  Air  Mass 


The  air  mass  element  models  the  momentum  flow  of  air  through  a  helicopter  rotor 
disk.  For  this  element,  the  rotor  is  assumed  to  be  an  actuator  disk,  and  the  flow  field  a 
cylindrical  region  surrounding  the  disk  (fig.  15).  The  state  vector  for  the  air  mass  element 
is  made  up  of  the  generalized  coordinates  for  a  single  air  node. 


Figure  15.  Air  moss  element  flow  field. 


Steady-State.  Consider  the  air  flowing  ste  lily  through  a  rotor.  Reference  2?  shows 
that  the  thrust  dT  acting  on  a  differential  annurus  of  the  rotor  (fig.  16)  is  related  to  the 
induced  velocity  v  via  a  momentum  balance  such  that 

dT  —  Airparv  \V  +  v\  dr  (8.2-1) 

where  r  is  the  radial  coordinate  of  the  rotor  and  T  is  the  velocity  of  the  rotor  relative  to 
still  air  (V  is  positive  when  the  rotor  is  moving  in  the  positive  x j  direction).  The  use  of 
the  absolute  value  of  the  sum  of  the  velocities  V  +  v  assures  that  the  differential  thrust  dT 
has  the  proper  sign  under  all  operating  conditions.  Integrating,  the  total  rotor  thrust  is 

R 

T  =  Anpa  I  v|V  +  v|r<i7’  (8.2-2) 
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Figure  16.  Air  mass  element  differential  annulus. 


The  virtual  work  done  by  the  thrust  on  the  air  is 


fR 

8W  =  4rcpa  I  vSR  |V  +  v|  r  dr 


(8.2-3) 


where  8R  is  the  virtual  displacement  of  the  air.  The  right-hand  side  of  the  equation  for 
the  virtual  work  can  be  discretized  by  letting  v  =  Of  -f  7j j*r  and  8R  =  8Pf  +  r8<j)fr.  Then 

fR  . 

8W  =4tt pa  j  (Of  +  ifrr)  |  V  +  Of  -f  |  {f>Pf  +  rHw)r  dr 


tR 

:8Pf4*pa  J  (Of  +  jfrr)  \V  +  Of  +  7,V|  r  dr+ 
8(j)^r4-npa  j'  (Of  +  7iV)  \v  +  Uf  +  7irT*| 7,2  dr 


(8.2-4) 


Note  that  while  the  coefficient  of  8Pf  in  equation  (8.2-4)  is  equal  to  the  rotor  thrust,  the 
coefficient  of  8<j>fr  has  the  dimensions  of  moment  but  no  clear  physical  significance. 

The  contributions  to  the  6W  (applied  loads)  side  of  equation  (8.2  4)  are  determined 
from  blade  element  theory,  and  are  obtained  by  summing  the  contributions  from  each  of 
the  aeroelastic  beam  elements  that  make  up  the  rotor. 

Dynamic.  Simple  models  for  the  induced  inflow  dynamics,  such  as  the  one  introduced 
in  reference  23,  have  been  shown  to  improve  the  accuracy  of  mathematical  models  of  heli 
copter  rotor  dynamics.  The  velocity  of  the  air  mass  is  idealized  as  consisting  of  a  spatially 
and  temporally  uniform  freestream  velocity  V ,  which  is  augmented  within  a  cylindrical 
region  by  the  steady-state  inflow  Of  induced  by  the  rotor  steady-state  thrust,  and  by  the 
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infinitesimal  dynamic  perturbations  to  the  inflow  induced  by  dynamic  perturbations  to  the 
thrust,  roll  moment,  and  pitch  moment  of  the  rotor. 


For  a  differential  annulus  of  a  rotor  disk  through  which  air  is  flowing  unsteadily, 
the  momentum  balance  can  be  expressed  as  a  system  of  first-order,  integro-differcntial 
equations. 


6W  = 


2rr 

2 pav \V  -1-  v|  8P  dxp  dr- {- 
pav6P  dVe{{ 


(8.2-5) 


Vcfi 


In  order  to  intermix  the  air  mass  terms  with  the  structural  generalized  coordinates  in 
a  single  set  of  second-order  equations,  the  perturbed  air  mass  generalized  velocities  are 
expressed  as  the  time  derivatives  of  generalized  coordinates. 


v  =  UA  +  7 frr  4-  P j  -  <p12rsinxp  -f  <p13rcos  xp 


(8.2-6) 


•_  A 

where  P j  is  the  vertical  component  of  the  perturbation  of  the  induced  inflow  velocity 

'A 

component  at  the  center  of  flow,  <^12  and  <p]3  are  the  flow  gradients  at  the  center  of  flow 
in  the  x2  and  v3  directions,  respectively,  and  xp  is  the  azimuthal  coordinate  of  the  rotor, 
measured  as  a  right  handed  rotation  about  the  Xj  axis  from  the  X3  axis.  The  flow  direction 
is  assumed  to  be  positive  along  the  x  1  axis. 

In  addition,  ■*  irtual  displacement  of  the  air  inside  the  cylindrical  flow  field  is  assumed 
to  be 

8P  =  8PA  —  8<pf2r  sinxp  +  6(pA3r  cosxp  (8.2-7) 

where  8PA  is  the  vertical  virtual  displacement  of  the  air  at  the  center  of  flow,  and  8<pA2 
and  8(pA3  are  the  cyclic  virtual  displacement  components  at  the  center  of  flow. 

Now,  consider  the  expression  u|Ir  +  uj,  where  v  ~  v  +  v(t).  In  seeking  the  linearized 
perturbation  of  such  an  expression,  if  V  +  v  ~  0  then  vjT  +  v|  —  (u  +  v)|v|.  Since  there  is 
no  linear  contribution  in  this  expression,  it  may  be  assumed  to  be  zero.  Now,  define 


{-|-1  for  a  >  0 

0  for  a  =  0  (8.2-8) 

- 1  for  a  <  0 


Then, 


v  |F  -|-  v|  =  [(V  -!■  v)v  I  (V  -f  2v)v]sgn(V  -I-  n)  (8.2-9) 
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Since  only  the  linear  perturbation  dynamics  are  pertinent  to  this  problem,  the  contri¬ 
bution  of  the  change  in  momentum  per  unit  area  term  from  equation  (8.2-5)  is 


nR  f2ir 


■.A  -.A 


If  2poK^  +  2u)sgn(F  -f  v)(Fj  -  ^12rsin^  +  ^13rcos  il))(6P*+ 
6<f)A2rsimf>  +  6 $f3r cos  if))  drj}  dr  = 

J  par{V  +  2v)sgn(V  4-  v)[P1  SP *  +  +  ^13^3)]  dr 

•  A 

Pj  6Pf  47t pa  j  (V  -f  2v)sgn(F  +  v)r  dr+ 

^Pai^n^n  +  ^i3^n)  +  2v)sgn(Vr  +  v)r3  dr 


(8.2-10) 


The  contribution  of  the  volume  term  from  equation  (8.2-5)  is  the  virtual  mass- virtual 
inertia  effect  as  calculated  in  reference  23. 


8  P 


f(R3  -  e3)J>  V,-*  +  “(a5  -  +  ilM,) 


(8.2-11) 


From  this  development,  the  coefficient  matrix  for  the  generalized  accelerations  may 
be  defined  to  be 


M  = 


8p„P3 


!-(ir 

0 

0 


R») 


0 

0 


(8.2-12) 


and  the  coefficient  matrix  for  the  generalized  velocities  may  be  defined  to  be 


C  =  2n  pa 


2  JtR  9?  dr 
0 
0 


0  0 

JR  9r 3  dr  0 

0  jf1  gr 3  dr 


(8.2-13) 


where  g  =  (V  +  2u)sgn(Vr  4-  v )  and  v  =  U f  +  r^fr. 

To  eliminate  all  periodic  coefficients  in  the  equations  of  motion,  and  to  assure  the 
existence  of  a  steady-state  solution,  the  air  mass  element  degrees  of  freedom  must  be 
inertial.  In  addition,  the  flow  direction  must  be  coincident  with  the  steady  state  spin  axis 
of  the  rotor  and  the  gravity  vector,  if  gravity  is  included  in  the  model. 
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8.3.  Aeroelastic  Beam 


The  aeroelastic  beam  element  is  designed  to  model  a  beam  undergoing  small  strains 
and  large  rotations,  and  for  which  shear  deformation  and  warping  rigidity  may  be  ignored. 
A  model  of  this  type  is  developed  in  reference  24,  which  formulates  the  nonlinear  beam 
kinematics  and  applies  them  to  the  dynamic  analysis  of  a  pretwisted,  rotating  beam  el¬ 
ement.  The  kinematic  relations  that  describe  the  orientation  of  the  cross  section  during 
deformation  are  simplified  by  systematically  ignoring  the  extensional  strain  compared  to 
unity.  The  only  restriction  on  the  magnitudes  of  the  orientation  angles  used  in  describing 
the  cross  section  orientation  is  that  they  remain  less  than  90°  .  All  influences  of  warp  other 
than  warping  rigidity  are  retained.  The  beam  cross  section  is  not  allowed  to  deform  in  its 
own  plane.  The  static  equations  from  reference  24  are  used  without  simplification;  the  dy¬ 
namical  equations  are  linearized  relative  to  static  equilibrium.  One  noteworthy  feature  of 
the  derivation  of  the  equations  in  reference  24  is  that  the  common  practice  of  using  an  or¬ 
dering  scheme  has  been  abandoned.  Thus,  all  higher-order  terms  (within  the  assumptions 
above)  are  retained. 

In  the  following  sections,  the  details  of  the  derivation  of  the  equations  for  the  aeroelas¬ 
tic  beam  element  are  presented.  First,  a  synopsis  of  the  basis  under  which  the  governing 
equations  of  the  beam  are  derived  is  given.  Next,  the  equations  of  motion  for  the  beam 
element  are  derived  in  terms  of  the  frame,  air,  bending,  extension,  and  torsion  degrees  of 
freedom.  These  equations  include  contributions  from  beam  elasticity,  inertial  and  gravi¬ 
tational  forces,  and  aerodynamic  forces.  Then,  the  discretization  of  the  beam  degrees  of 
freedom  is  presented  to  show  how  the  beam  displacements  are  transformed  into  the  beam 
generalized  coordinates.  The  final  two  sections  describe  the  transformation  from  root  and 
tip  node  degrees  of  freedom  to  beam  generalized  coordinates,  and  the  transformation  from 
beam  generalized  forces  to  root  and  tip  node  forces  and  moments. 

8.3.1.  Basis  of  the  Governing  Equations 

Consider  the  beam  clement  shown  in  figure  17.  The  element  frame  is  denoted  by 
F,  and  the  root  and  tip  nodes  are  denoted  by  R  and  T,  respectively.  The  addition  of 
primes  and  double- primes  signifies  the  static  and  perturbed  dynamic  states,  respectively. 
It  should  be  noted  that  F"  and  R  are  coincident  with  each  other  and  *hat  their  coordinates 
line  up  with  the  principal  axes  of  the  root  end  of  the  undeformed  beam  element  with  the 

undeformed  beam  lying  along  b3  .  Similarly,  T  is  at  the  tip  of  the  .mdeformed  beam 
element  and  its  coordinate  directions  lie  along  the  principal  axes  for  the  tip  cross  section. 
The  air  node,  denoted  by  A,  must  be  included  in  the  problem  so  that  the  influence  of 
aerodynamic  forces  on  the  air  node  generalized  forces  can  be  determined  and  so  that  the 
influence  of  perturbations  of  the  air  node  generalized  coordinates  can  be  determined  for  the 
generalized  coordinates  of  both  the  beam  and  the  air  node.  The  position  and  orientation 
of  A  are  inertially  fixed. 
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Figure  17.  Aeroelastic  beam  element  (undeformed  with  pretwist). 


Interior  displacements  of  the  beam  are  represented  by  four  functions  of  the  axial 
coordinate  xy.  ui  and  0 3.  Bending  is  described  by  and  a  2 ,  axial  displacement  by  113, 
and  torsion  by  63 .  These  functions  are  discretized  in  terms  of  standard  cubic  and  linear 
polynomials  so  that  the  generalized  coordinates  at  the  root  and  tip  of  the  beam  can  be 
related  to  the  nodal  displacements  and  rotations.  In  addition,  however,  there  are  also 
generalized  coordinates,  called  internal  degrees  of  freedom,  associated  with  higher  order 
polynomials. 


8.3.2.  Beam  Elasticity 

The  derivation  of  the  equations  to  calculate  contributions  of  the  elastic  deformations 
of  a  straight,  pretwisted  beam  follows  the  derivation  presented  in  reference  24. 

Steady-State.  The  elastic  beam  equations  for  a  beam  in  equilibrium  are  derived  from 
the  variation  of  the  strain  energy 


SU  =  /  {GtZaStza  +  226336633)  d®3 
</o 


(8.3.2-1) 


where 

<=3i  =  (Ai  -  £2)  («3  -  Q1) 


^32  =  (^2  +  £1)  («3  -  Q') 

633  =£33  +  -  £l^2  +  -  (£j2  +  ^22)  (^3  -  0')2  +  (&Ai  -  £1 A2 )  («3  —  O')  O' 
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where  9(x3)  is  the  pretwist  angle  (fig.  17),  with  9(0)  =  0,  and  ()'  =  <£(  )/dx 3.  The 
generalized  strains  are 

633  —s'  —  1 

s'2  =u'2  +  v!2  +  (1  +  U3)2 

Ki  =(Ci2tt"  -  Cnu2)/C33  (8. 3.2-3) 

K2  =(C22u'(  -  C21u'l)/C33 

ki  =«;  -  +  CsiM")  /Cjj 

where  C  =  CP  F' ,  the  direction  cosines  of  local  principal  axes  relative  to  the  static  frame 
orientation.  The  elements  of  C  may  be  expressed  in  terms  of  Tait-Bryan  orientation  angles 
(orientation  angles  of  type  body-three:  1-2-3)  as 

Cll  =C2C3 

C\2  —s2Ci  +  Sia2c3  (8.3.2-4a) 

C13  =53^1  -  CiS2C3 

C21  =  -  C2$3 
^22  =C3Ci  —  S1S2.S3 
C23  =C.3Sj  +  CjS253 

(8.3.2-4b) 

C31  =Uj 
C32  =U[ ; 

C33  =(l-u'2  -u'22)’ 
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where 


*1  =  -  t4(i  --  u'i2)~^ 


S-i  =Uj 

a  3  =  sin  0j 

(8.3.2-S) 

ci  =(1-V)* 

Cj  =(1  —  323)  * 

C3  =  COS  03 

After  integrating  over  the  cross-sectional  area,  the  variation  of  the  strain  energy  is 
obtained  in  terms  of  the  stress  resultants  JF3,  Mi,  M2,  and  Mj. 


=  I*  (Fits' 
Jo 


+  M\Sk\  +  M28K2  +  MjSkj)  dx 3 


(8. 3.2-6) 


where 


F3  =EoC33  +  E2K1  —  E\K2  +  jt32  -j-  Dq6't3 
Mi  —E2Z33  ~f*  -Tj/c  1  H — - 1-  D20't3 

It 


D  _  2 

M2  =  —  EiC 33  +  I2K2 - ~Z - D\9't3 

/t 


(  B3T32  3ZM' 

M3  =  {j  +  lihs  +  Kj  ~  E1K2  -I - — - 1 - — — 


-r3  +  D46'2^J 


(Dq€33  +  D2K 1  —  DiK2)6' 


(8.n  ::~7) 
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where  r3  =  n3  —  6'  and  the  section  integrals  are  defined  as 

E0  =  JjEdA 

Ei=JJ  E^dA 

E2  =  JJ  E(2dA 
h=JJ  E(22dA 
h=  JJ  EtfdA 

h=h+h 


J  =  J J  G[(Ai  -  fo)2  +  (A2  +  (i)2)dA 
B1=  J J  JBfctfi*  +  (22)dA 
B2  =  JJ  +  6*)«m 
B3  =  JJ  Eitf  +  tffdA 

Do=JJ  £(fcA,-*iA2)cLA 

Di  =  J J  E(i(£  2A1  -  (2X2)dA 

D2  =  J J  E£2(£2\i  -  £i\2)dA 

D3  =  IJ  JB(6*  +  6a)(6A,-e,A,)iU 
DA  =  j  j  E{(,2\i-ix\2)2dA 


(8.3.2-7b) 
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Here,  Eu  is  the  axial  rigidity;  E\  and  E2  are  the  first  flexual  moments  about  the  local 
£1  and  £2  axes,  respectively;  I\  and  I?  are  second  moments  ( bending  rigidities)  about  the 
local  £1  and  (2  axes,  respectively;  and  J  is  the  Saint-Venant  torsional  rigidity. 

The  variations  of  generalized  strains  can  be  expressed  in  terms  of  the  fundamental 
variables  as 

6s'  -~-6ui 

OUi 

(8.3.2-8) 

~~Q^n  ■*"  ^3*^3  +  ^rK  +  c3«a*a^3 

and  the  variation  of  strain  energy  as 

w  “  jf  [(*&  ‘ +  M‘S)  ^  +  + 


8k- 

Mig^Sul  +  mSO'A  dx 3 


where 


ds'  63i  +  u'j 
du’-  s' 


du1;  C33 


dny  C22 

dv.”  C33 


(8.3.2-10a) 


9k 3  _  O33 2 

9u'{  C33(l  -  C3i2) 


c332  V  1  - 


C22C32  ,  CJ2C3A  UjCiiCsiCf, 


C312  ^ 


^333(1  -  C?312) 
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^*2  _  ul  (  CliCM  ^22^31  \  _  U;|C2lC3iC322 
^ui  C33 2  \1  —  Cji2  G33  /  C'333(l  —  C3i2) 


9K3 

du'j 


-C31C32U" 
Cf333(l  -  C3i2)2 


(2C332  +  C312  -  Csi4)  - 


«a(l  ~  <?32 2 ) 
C333 


dn  i 
du'2 


-Cj  1 
C333 


(«jC31  +  U2C32) 


(8.3.2-10b) 


Dynamic.  Since  the  explicit,  analytical  derivation  of  the  elastic  stiffness  matrix 
would  be  exceptionally  tedious  and  lengthy,  GRASP  generates  it  numerically.  This  is 
accomplished  by  taking  the  Jacobian  of  the  function  that  calculates  the  steady-state  elastic 
loads.  Because  of  the  necessity  of  calculating  an  accurate  stiffness  matrix,  the  algorithm 
used  to  calculate  the  Jacobian  uses  n  two-point  central  difference  scheme  plus  a  generalized 
formulation  of  Richardson  extrapolation. 


8.3.3.  Beam  Inertial  and  Gravitational  Forces 

The  generalized  forces  resulting  from  motion  of  the  aeroelastic  beam  relative  to  an 
inertial  frame  are  also  determined  following  reference  24.  Warping  dynamics  are  again 
ignored.  The  derivation  is  based  on  the  work  done  by  inertial  and  gravitational  forces 
moving  through  a  virtual  displacement.  The  work  is  calculated  by  taking  the  scalar  product 
of  the  gravity  minus  the  acceleration  of  a  generic  point  P  (fig.  17)in  the  beam  interior 
(ref.  24,  eq.  32),  with  the  virtual  displacement  of  the  same  point  (ref.  24,  eq.  34),  then 
integrating  the  result  over  the  beam  length. 

Steady-State.  For  a  beam  element  in  equilibrium,  the  virtual  work  is 
£W  =  /  (6utFf',  +  Sif)piMp,)dx3+ 

J  0 


(8.3. 3-1) 


where 


=m(gF,  -  <')  -  hFF:,nF;,'R’?,r‘ 


and 


Mpi  j  =T7ljCrpfS  — *lfip»2^P'3 
Mpi  2  —  —  T7liGp»3  4"  *2^jP'i^P'3 


A/pij  =m\Gp'2  —  m^Gpii  +  (»i  —  t2)^p'i^p<2 

iZjsvf  =m(x3^3i  +  n»)  4-  miCu  4-  miCn 

r 

c?p*  —qf>  —  Api^  — 

u  =Rp'lp 


6u  =6Rp,p 


6lppi{  —SjiSOj  +  g^&u'a  =  typj 

Supi  =SRp',F 


6^f>  —Sipp, 


p'p 


The  section  integrals  are 


m.\ 


7712  = 


*1  =  j j  P'tldA 
*a  =  j  J  P&dA 


(8.3.S-2) 


(8.3.3-3) 


(8. 3.3-4) 
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Dynamic.  For  a  linearized  perturbation  about  the  equilibrium  solution,  it  is  possible 
to  express  the  equations  of  motion  in  a  matrix  format  such  that  the  virtual  work  per  unit 
beam  length  is  given  by 


f  SuF"i  ^  1  UF"t  UF"i  fUF"* 

S^F"i  I  0F"t  0F"i  &F"i 

-6W  =  <  6ui(x 3)  J  <  [M]<  £<(*3)  ►  +  [C] <  ^(*3)  >  +  [Jf]<  «.(« 3)  >  -{<?}> 

^ua(*s)  **{**)  ^(*3)  “a(*3) 

(8.3.3-5) 

where  the  components  of  the  generalized  force  vector  Q  are  the  same  as  the  static  general¬ 
ized  forces  (see  the  previous  section)  and  the  coefficient  matrices  M,  C,  and  K  are  defined 
on  the  following  pages.  The  M  coefficient  matrix  is  defined  as 


SuF«i  row: 


UF"j  column:  mSij 


9 F" j  column:  eijk(mRp,F  +miC(jt  F  ) 


Uj  column:  mSij 


(8.3.3-6a) 


v-'p  column:  efczmmmC£ 


p>  pi  0k\ 


9 3  column:  e*3imiC£  F 


95 


Stf>F"i  row: 


«f»*j  column:  —  ^(mJCfn  +  mjC^  ) 


tP'F' 


Bp"j  column:  6ij(mRpi F  Rp'k  +2 to/C|*  F  )— 


*P'F’  nP'F1  • 


_  nP’F'  nP'F'  m  ffiP'F'  nP'F'  ,  >oP'F'  nP'F'' 
mRpn  Rpij  —  mi{Cti  lipij  +Otj-  iijwi  , 


itj  column:  ~  <ijfc(m^F'fc  +  ml^ik  F  ) 


•j'  ,  _  ( ,-,P' F'  fiP' F'  siP'F'rP’F's  nP'F'  ■ 

Up  column:  Cim  -  Cfem  CH  )*iF'm  ^-+ 

•  rp'F'^!±  L,„rp'F'—  • - 

llCl<  2  2<  au^+t3C3‘ 


column:  mt(C7^  F  CFt  F  -CFtF  Cj£  F  )i?£.f  + 


Sui  row: 


ur"j  column:  m0y 


tP'F' 


0F"j  column:  eijk{mRp,F  +  miCfk  F  ) 


Uj  column:  mtfjj 


tip  column: 


5/Ci 

0U" 


03  column: 


tk3imC&F' 


(8.3.3-6b) 


(8.3.3-6c) 
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fiu'a  row: 


pi  pi  Ski 


uFnj  column:  CkimmmCfc  — 


h"j  column:  mt(c£F'  CfmF’  -  Cf-F<  C£^'  )Rp,^  ~-4 

•  &K1  siP'F'  ,  •  siP'F'  ,  :  &k3  fiP'F' 

"wc'i  +”wcv  +,Ia<05' 


pi  pi  8k\ 


uj  column:  €kimmmCkj  — 


.  8k\  8ki  .  &K2  8k2  .  8k3  8k3 

Up  column:  t,  ^  ^  +  *2  g~^  ^  +  *3  gu„  du„ 


03  column:  13 


8*3 

d< 


03  row: 


ti.F"j  column:  Ck3imiCFjF 


«r»,  column:  mt«F’c£'F'  -  C^'C^)R^  +  i5Cg 


tP'F' 


iij  column:  tkumiCkj 


pi  pi 


-it  ,  .  d/c3 

column:  13^77 

OUp 


63  column:  13 


(8.3.3-6d) 


(8.3.3-6c) 
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the  C  coefficient  matrix  is  defined  as 


SuF"i  row: 


•  S’*  T 

UF"j  column:  —  2tijkm\lp,k 

&F"j  column:  26ijflpij[{mRFiF  +  miCjitF  )— 
2fi£!*(m.RF!f,  +  mkCu  F  ) 

•  p*  r 

uj  column:  -2 eijkm,Q,F'k 


l>  1  O  nF'IfnP'F'  nP'F1  r>P'  F'  r*P'  F'  \&Km 

Ug  column:  2mi>n F'li^mi  ^ki  . ki  ^ ml  ) 


H 


6tj)F"i  row: 


k  column:  2m*nF:/(CF/F'CF/F^  -<f'c£f') 

up"j  column:  —  26ijflFt£(mRpiF  +  F  )+ 

2n£!-(mi?£!F'  +  mjfeC'£'F') 


i  i  ciFf  I  jdP9  F1  jyP1  F9 

9F"j  column:  —  2e;jfcT7iHpij  Rp>k  F'l 

2««.mmn?:/(cr,F,ii?:r + czr^r')- 

2eijkHQ’F'l 

■iij  column:  —  2SijUFik{‘fn^F,k  +wii Q*  ^  )"b 
2n^!F(7ni2^! j5’*  +  tt ikCFjF ) 


up  column:  2efc(mmm(ftF,fCF„F  -nF»f,C'jf<F  )Rp’n 

P'F'\r,P'F'  ®Kl 


2flF,^»(efc/2tlC^nF  +  tkll^ClnF  )Cki 


du"p 


03  column:  2efci3m{(fi£,  j CFmF  -  Up'L^ki F  )RF'm  + 


(8.3.3-7a) 


(8.3.3-7b) 
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Sui  row: 


it pn  j  column:  -  2c ijkmUp,{ 

Bpnj  column:  2 6ijSlpi[{mRp,F  +mj Cfc  F  )— 
2Slp,j(mRp,f'  +micCFlF) 


F'l 


iij  column:  —  2 


column:  2 -  C^'C^')^ 


9,  column:  C^'  -  C^'cD 


Su'  row: 


column:  -  2m»n£f(C#-’Cf('F'  -  CfcJl'')8"' 


F  r<P  F  \__HL 

dua 


\F'Ir>P'F'  nF'Ir'P'F'\nP'F'  &Kl 


eF»j  column:  -  2c*immm(fl£,JCfcnF  -  nF>nCkjF  )Rp,F  — H 


2nf,U^cC  +  'mcCK*' 


iF'r(r>P'F'  nP'F'  nP'  F'  nP'  F1  \&Km 

*•*'  3u" 

a 


tij  column:  -  2mhn£,j(CFf  Ckl  1  -Ckj  Cml  ) 


Uq  column:  -  2fi^^(efcl2iiC2p^F'  +  iktli2C^F  )—• 


dua  dup 


6,  column: 


(8.3.3-7c) 


(8.3.3-7d) 
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66 3  row: 


UF»j  column:  —  2mkUF,f(CPjF  C pt  F  —  C^F  CPjF  ) 


$F"j  column:  —  2e*3imi(ft£,^CFmF  —  ^F'm^kj1^  )^F'm~ 


iF'f  ^P'F'mjP’F’ 
'm 


2nFF:i(ilcFiF,cF;F'  -  i2crk*c2y) 

*j  column:  -  2 CP'F'  -  C^'C?/') 


P'F'nP'F •' 


(8.3.3-7e) 


*)j  column:  -  2n£i(;.c£F'f|  -  <.05''f|) 


6*3  column:  0 

and  the  li  coefficient  matrix  is  defined  as 


foiF»  row: 


up"  column:  mUp^Up,1 


Op"  column:  -  gp<  +  ( Up^UpjR P,F  ) 


-  Up^Clp/rhF' 


nF'lnF'j- 


-nFF:inFp:iRpF:F'}+ 

(8.3.3- 


u  column:  mil p,1  Up,1 


6pi  column:  —  Up^Up^mp'O1^  p 
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Sippn  row; 


uf"  column: 


-  m(i£'  -  -gF.  +  r$,f‘  )■  -  - 

(0F!J  flF^mp')  +  ThFihF,ItlF,1 


Of"  column:  7n{i?p/F  {Apt1  —  5F')+ 

R^'i^n^R^')'  -  nFF:inFF:iRpF:F')}+ 

mHip  -  gF>  +  {fl^n^Rpf'Y- 

n£TnF:rRpF:F')+ 

RF'F'  [(ngl'n^WO"  -  n^/nf!/mF']+ 

cF,p,(#p,f2g/-^:jtp.^:/)cp'F' 


(8.3.3-8b) 


u  column: 


-  m[Ip  -  gF.  +  (f 

(iip^iip^mF')  +mF»nF,/nF»/ 


dpi  column:  -  i2piF  fiF/flF.JmF*CF  P  — 
(rhFinFF:ihFF:iRpF:F')~cF,p,+ 
mFi(nF:inF:iRp:F,)~cF'p'- 
(mF^Ap!1  -  5F')1  CF  P  +  mp^Ap,1  -  gF>)CF  p  -f 

cF'p'(ifp.n£/  -  nflVn?;1)7, 

5u  column: 

up"  column:  mOp/Hp/ 


dF"  column: 


m[AFF,1  -  gF'  +  (n pu^R^')'  -  ftpn 

(np/np/mp.)  -  nFiJ j2pl/mp« 


F'lffP'F'  i , 
pi  xtjri  I  -r 


(8.3.3-8c) 


u  column: 


mUpJflp,1 


6pi  column:  —  fip/Hp^mp*^  p 
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Sifrpi  row: 


uf"  column:  Cp'F>  mFittp'^^F'1 


where 


6p»  column:  Cp  F  mp^Api1  -  gF')  +  Cp  F  thp' [({Ip'1  (ipi1  Rf>F  )  - 

,F'/oF'/6P'F'i  ,  _  nF'/ •  . oF'/^P’F' 


nynyRfr*  }  +  {HP.ny  -  nyip-nj/jc1 


(8.3.3-8d) 


u  column:  Cp  F  mpittp^ilp,1 

8p>  column:  Cp  F' rhpilAp,1  -  gp>  +  (flFiIflF>IRpiF')  ]CFP  + 


Hp.nfr1  - 

oF7* 
*lpi  ip1 

I 

HP , 

=i  p*  rip 

/ 

u 

0 

o' 

iPi 

= 

0 

*2 

0 

0 

0 

*3 

(8. 3. 3-9) 


In  the  foregoing  matrices,  m  is  the  running  mass  per  unit  length,  and  ma  is  the  first 
mass  moment  about  the  £a  axis.  The  last  block  row  associated  with  8^>p,  p  is  used  to 
obtain  the  terms  associated  with  8u'a  and  80 3  by  substitution  from  the  equations 


8tJj 


p"p'  _ 

P'i  - 


htfQj  +  {-7rh  + 


d2K{ 

du'^du'p 


d2Ki 

fTMdk 


U*  + 


83)8u'a 


(8.3.3-10) 


and 


nP"P'  r  0  ,  ®K'  ~.l 

0p>i  -  O^Oi  +  fanua 


(8.3.3-11) 


The  8rj)P<  row  block  matrices  must  then  be  pre-multiplied  by  RT  and  the  6P>  column  block 
matrices  must  be  post-multiplied  by  R,  where 


Ria  = 


8k  i 

OK 


(8.3.3-12) 


Ri3  =8ii 
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The  geometric  stiffness  matrix  KG  is  added  to  the  6if)p>  rows  and  8p>  columns  where 


vG  - 
■h-Q/3  - 


82k  i 
du'^du'p 


*  £  =  -  Afft 


d2Kj 

du'^dOs 


(8.3.3-13) 


The  matrix  KG  comes  from  the  last  two  terms  in  equation  (8.3.3-10),  which  are  commonly 
called  the  geometric  stiffness  terms. 


8.3.4.  Aerodynamic  Forces 


The  aerodynamic  forces  acting  on  the  aeroelastic  beam  element  are  determined  from 
a  quasi-steady  adaptation  of  Greenberg’s  thin-airfoil  theory  (ref.  29).  Before  the  theory  is 
discussed  in  detail,  two  new  sets  of  axes  must  be  introduced  for  the  purposes  of  defining 
the  directions  in  which  the  lift  and  drag  forces  and  the  pitching  moment  act.  In  figure  18, 

the  Z  axes  are  associated  with  the  zero-lift  line  for  the  airfoil  section  with  the  vector  b2 

-  z 

along  the  zero-lift  line  toward  the  trailing  edge.  The  vector  b3  is  along  the  beam  axis  but 
in  a  direction  such  that  a  dextral  rotation  of  the  airfoil  section  about  ^iis  vector  results 
in  an  increase  in  the  angle  of  attack.  Then,  being  a  dextral  system,  turns  out  to  be 
normal  to  the  zero-lift  line  (and  nominally  in  the  direction  of  positive  lift  for  the  section). 


-w 

The  other  set  of  axes  is  the  so-called  wind  axes  W.  For  these  axes  the  base  vector  fe3 
is  identical  to  b3  .  The  base  vector  b2  is  located  along  the  relative  wind  vector  (in  the 

-IV 

direction  of  drag)  and  6j  is  in  the  direction  of  lift. 


Zero-lift  line 
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The  Z  basis  and  the  P  (principal  axes)  basis  convect  with  the  blade  cross  section,  and 
are  related  by  the  direction  cosine  array  CZP  —  Cz  p'  =  Cz  p  " . 


Czp  = 


cos  8X  sin  0X  0 
-asin8x  crc.o$6x  0 
0  0  a- 


(8.3.4-1) 


where  a  =  +1  if  a  dextral  rotation  about  63  results  in  an  increase  in  the  angle  of  attack 
and  cr  =  -1  if  a  dextral  rotation  about  63  results  in  a  decrease  in  the  angle  of  attack. 
The  wind  basis  W  is  related  to  the  Z  basis  by 


cos  a  —  sin  a  0 

Cwz  =  sina  cos  a  0  (8.3.4-2) 

°  0  1 

where  a  is  the  angle  of  attack.  Then,  Cwp  =  CWZCZP . 

Point  Q  is  the  quarter-chord  point  of  the  cross  section,  about  which  the  aerodynamic 

forces  and  pitching  moment  are  calculated.  The  offset  position  of  Q  relative  to  the  origin 

+ z 

of  the  local  principal  axes  P  is  R^P b2  . 

Consider  the  wind  velocity  vector  at  the  perturbed  position  of  the  aerodynamic  center 
Q".  WQ"  is  calculated  by  subtracting  the  inertial  structural  velocity  at  Q"  T)  from 
the  inertial  air  velocity  at  Q"  (U®  f),  where 


=  -(Iff  +  rf*  +p*  +  n%}£,  + 


(8.3.4-3) 


V f  =(ftF'J  +  ifF'  +iLPF)RQP  +  (nF>I  +  qf"f)rpf"+ 


jF"  F'  ,  F'bF"F'  ,  JbF>I 


R  +  fi  R*  *  +  *  R  +  *R 


The  relative  wind  velocity  components  in  the  Z"  basis  are  then 


wgli  =  -  C?“A(VA  +  rjA  +  P*  +  R%Ain  +  R%Ainh 

((c2"F"o?:')'  +  (cz"F'V.f  +  {c^<-'Cp')')rY- 
((c2"fx:')' + (cz"F"'eF„)')cz”,”'Rp;:r- 
cz"F"  i  -  i(c2"F’’uf,I)'}cz“F"iF„- 
(cz"F"dF,,)'cz"F"uF„  -  cz"F"iF..  -  cz"F"v£I 


(8. 3.4-4) 
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end  the  local  air  flow  velocity  gradient  is 


aw$, 

z"‘ 2  9i?ff 


=  -  cfA?£«2Lt£»L*2LfA  +  CA*£  +  CA*tft)+  <8-3-4-5) 


«  +  cf/'i;  +  c§< 


where 


If  the  time  derivatives  in  equations  (8. 3.4-4)  and  (8.3.4-5)  are  replaced  with  variations, 
the  relative  virtual  displacements  and  rotations  of  an  element  of  air  with  respect  to  the 
structure  are  obtained. 

*S%: ,  =  -  Cf,"A(SP,A  +  R%A6tf2  +  R%aW?3)- 

[(cz"F"s^,y + (e2"p"^?::p,)-]jlB«'>- 
(c^'H^iiCf^R^r  - 

Cjf  p"« uF„i  -  (Cz"f"6^f.,)C2"f"Rf;„f‘ 

n%:,  =  -  cx*cg*s+t,  -c?;ac%ah  a3 + c?;,F"6^i+cf;,F"6tf:F' 

(8.3.4-6) 

The  relative  wind  velocity  magnitude  and  components  are  time- dependent  quantities. 
For  the  magnitude  note  that 


for  which  the  static  part  is 


and  the  dynamic  part  is 


w! =(w$:j+(wg:1f 
w*  =  +(w$ty 


w- 


w<?'  t yQ"  i  urQ'  i,M" 

* '  Z'l  ' '  Z"  1  +  S' 2U  Z" 2 

W 


(8. 3.4-7) 
(8. 3.4-8) 
(S.3.4-9) 


Likewise,  the  angle  of  attack  is  a  time- dependent  quantity.  In  the  equations  written 
below,  it  is  necessary  only  to  develop  the  static  part  and  the  linearized  dynamic  perturba¬ 
tion  part.  Those  quantities  are  easily  determined  from  the  definition  of  a. 


tana  = 


ty*?" 

yYZ"l 

WQ" 

n  Zn  2 


(8,3.4-10) 


105 


The  static  part  is  simply 


i irQ' 

yyZ<  2 

(8.3.4-11) 

while  the  dynamic  part  is 

W^'  w Q"  —  W^'  IV®" 

vv Z'2n  Z"\  VYZ'lvyZnl 

flrJ 

(8.3.4-12) 

The  applied  force  is  assumed  to  be 

-IV"  -IV"  .7I» 

F  =  Ccb,  +Vb2  +Cncbf 

(8.3.4-13) 

and  the  applied  moment  is 

M  -  Mb  f 

(8.3.4-14) 

The  equations  governing  the  aerodynamic  force  components  are 

cc  =\p«w*cci  +  ^.c2waf„n 


V  =-p,W*cc4 

M  =i p,W*chm  -  ^,.c\WG%.n  +  wg:t  +  jG%:u) 
Cm  =  j/vW,',  +  C-G%,n) 


(8.3.4-15) 


Now,  all  of  the  quantities  that  are  needed  to  define  the  virtual  work  are  available. 

8W  ~  [  {-6S%'liFz"i  +  6T%',3M)dx3  (8.3.4-16) 

Jo 

Steady-State.  The  static  generalized  forces  can  be  removed  from  the  expression  for 
the  virtual  work  and  written  in  the  form  6W  =  fj  6qTQdx3,  where  5q  is 


0(P  12 

8q  =  8up"i 
6l{>F"i 
6ui 

Su'a 

‘60  3 


(8.3.4-17) 
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and  the  elements  of  Q  are 


8P*  : 

CeCl"" +VCI"" 

Hi- 

r(CeC"r'+VC *”•)  MC“ 

Hi- 

0 

Hi- 

0 

Supn  : 

CcC*fF'  +vc£,F' 

SrJ)F'i • 

MCfr+iRj'r'F^ 

Sui  : 

Ff 

8u'a  : 

c”(M-n%r*)%L 

80 3  : 

ClP(M  -  R%F*) 

(8.3.4-18) 
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Dynamic.  After  removing  the  steady-state  contribution  to  the  virtual  work,  the 
virtual  work  per  unit  of  beam  element  length  done  by  the  aerodynamic  forces  and  pitching 
moment  can  be  put  into  the  following  form: 


'8P* 

Htr 


i 


-sw  =  l 


Hti 
SuF"i  } 
Sui 


[A){ 


Cc 

V 

£>nc 

M 


\s< 

U03  f 


(8.3.4-19) 


This  equation  can  then  be  rewritten  in  terms  of  aerodynamic 


i  j  ,  auu  li.  ijuatuico, 


M  -AFH 


C  -AEH  +  AFG  +  BH  (8.3.4-20) 

K  =AEG  +BG  +  D 
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The  elements  of  A  are 


SPA  row: 

Cc  column:  —  A 

V  column:  -  G2['A 
Cnc  column:  —  C21 A 
M  column:  0 

6<f>Ar  row: 

Ce  column:  0 

V  column:  0 
Cnc  column:  0 
M  column:  0 

8(f)A2  row: 

Cc  column:  —  AR%A 

V  column:  -  C™  ARa2A 
Cnc  column:  -  CfiAR%2A 
M  column:  Cf  \  AC22  A 

5<j>A3  row: 

Cc  column:  -  AR%A 

V  column:  -  C2\  AR%A 
tnc  column:  -  Cfj  AR<\A 
M  column:  Cfj  AC22A 


(8.3,4-21a) 


(8.3.4-21b) 


(8.3.4-21c) 


(8.3.4-21(1) 
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8uF"i  row: 


Cc  column:  —  Cjf  F 

V  column:  —  Cjf  F 

Cnc  column:  -  F 

M  column:  0 

6if>F"i  row: 

j Cc  column:  tijkCyj  F  R%k 

V  column:  e.jfeC^y  F 

w  t pi  qJ 

Cnc  column:  e ijkCfj  lip, k 
M  column:  —  F 

6ui  row: 

tc  column:  -  C™  F 
V  column:  —  C™  F 
Cnc  column:  —  C ft 


(8.3.4-21e) 


(8.3.4-21f) 


(8.3.4-21g) 


(8,3.4-21h) 
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69 3  row: 


4  column:  C*'*  C*,PR% 
V  column:  Gjjf  z  Czp  R^p 


£nc  colunm:  C3P  R^P 


M  column:  —  C*p 


'33 


The  elements  of  B  are 


6Pf  row: 


tv'*' 


•noir^^+p^-) 

|tv?|2 


column: 


IV§ j  column 


Itvll2 


ruw: 


column:  0 


W^i  column:  0 


W§. 2  column:  0 


G§j2  column:  0 


S(j>i2  row: 


T'F^j  column: 


i^?p 


column: 


urQ  nZ'A  ftQ'A(r  2£ii 
UA2  (Lc  d& 


it”*' 


ult'Z’ 


(8.3.4~21i) 


(8.3.4-22a) 


(8.3.4~22b) 


(8.3.4-22c) 


G^jj  column: 


0 
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Wii i  row: 


W?,  column: 


!<l2 


\\'z2  column: 


G§  12  column:  0 


W^C^R^iCc^^  +  V1- 


(8.3.4-22d) 


\W$\> 


Supi'i  row: 


Wzi  column: 


Wz2  column: 


=V^(CJ^-±D?C« 

}W§? 


K2,C*A(C'2fc-  +  v- 


(8.3.4-22e) 


\W$\2 


Ggi2  column:  0 


SiJ)F»i  row: 


column: 


W%  column: 


in'll1 


in'll" 


i  ariW'i'  nr.w'x' 

t  r  ‘'O.f  ,  ir\OL>,i 


(8.3.4-22f) 


c  8  at 


G§j2  column:  0 


Sui  row: 


in'll* 


W%{  column: 


Wzt  column 


C?z  12  column:  0 


wS  C?'F'(C  ) 


(8.3.4-22g) 


in'll2 
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6u'  row: 


tlrQ  CZPH9P  8k *  (C  *  4.  7}ec»  *  •> 

*VZ2°33  ■*t^2  a*  r  u  §a / 

V|  p 


fWS9 


W§i  column: 


W%2  column: 


.virQ  C^p  R®p  C.  1  |  T>eg»  2  s 

nz 2  ou da  +  ^ aa  1 


|w?p 


<jfla  column:  0 


#03  row: 


w'g;c3Y-Ri2r(£e^ri-  i 


-.»v'z' 


W%t  column 


in'll2 


in'll2 


VTj-j  column: 


Gj„  column:  0 


(8.3.4-22H) 


(8.3.4-22i) 


113 


The  elements  of  D  are 
SPA  row: 

PA  column:  0 
<j>Ar  column:  0 
<j>A2  column:  0 
<j>A3  column:  0 


upiij  column:  0 


&F"j  column:  €kij(CcC™'F  +  VC™  F  )CF1A 
uj  column:  0 

u'0  column:  -  (klm{C.cC™‘ P>  +  VC^'P')C^A 
03  column:  -  ekn(CcC}l'p  +  VC3l  p  )Cpl  A 


d*m 


(8.3.4-23a) 
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5(fiiv  row: 


P/1  column:  0 
04r  column:  0 
<j>M  column:  0 
column:  0 
up-d  column:  0 
Opuj  column:  0 
uj  column:  0 
u'p  column:  0 
93  column:  0 


(8.3.4~23b) 
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8<j)y 2  row: 


PA  column:  0 
<j)Ar  column:  0 
<f)f2  column:  0 
column:  0 


uFnj  column:  -  (Cc C™'A  +  VC™'A)C£ 
h;,  column:  -  +  VC** 

c*rc[;*R%A)+ 


iW'A\r>F'A 


Uj  column:  -  {CcC A  +  VC™  A)C 


n 


4  column:  Cum[(CcC^‘  +  VC^Cg- ■ 

*nc$rct;*cg*  -  c$*cC*qni%± 


03  column:  ckl3[{CeC%' z'  +  VC%'z'){Cg- AC^AC^PRQZP- 

CfLPCr’AR%A)+ 

M{c?kpcr;Ac*A  -  c?;ac&ac2y)) 


(8.3.4-23c) 
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Hu  row: 


PA  column:  0 
<£4r  column:  0 
<j>A2  column:  0 
ififz  column:  0 


upnj  column:  —  {CcC^  A  H-  VC 2\  A)C*. j 


t  w'A\riF'A 


lr»i  column:  -  eiU\(CcCZZ'  +  VC^' )(cZ,*c£AR%Fk' 

cSTci i'^-V 

KAtrZ'F'  rF'ArZ'A  rZ'A^F'ArZ'F'w 

>/VH Cli  °Jfel  °23  “  ‘■'ll  °i3  C2fc  ;J 


iij  column:  -  (£cC1u/  4  +  A)Cf3 
:ijj  column:  cnm((£cC""2'  +  »C  z’  ^ c£ J cff  R^'- 

CgP<AR%A)+ 

M(czkpcg'Acz;A  -  Cf/Cf/Cf,")]: 


iW'A\r<F'A 


iZPriP'AnZ'A  n Z' A^P' A/nZPy] 

du* 


03  column:  c^-AT*'  +  VC%'*')(Cg- AC[,AC^PR^~ 

cHpcH,arqAa)+ 

M{c?kpcr;Ac*A  -  c?;ac&ac2y)) 


(8.3.4-23d) 
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8uF«i  row: 


Pj4  column:  0 

<frfr  column:  0 

<j>y 2  column:  0 

<j>}3  column:  0 

■UF'ij  column:  0 
Op"}  column:  0 
iij  column:  0 

u'„  column:  -  eUm{CcC] p*  +  VC%P'  )C,P  ^ 

83  column:  -  e*t3(£cCi*  p  +  VC^\  p  )Cfi  F 


(8.3.4-23e) 
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8if)F»i  row: 


.P/4  column:  0 
<f>ir  column:  0 
<f}^2  column:  0 
0^3  column:  0 


up"]  column:  tmki{CcC^mz  +VC2mz  )Cjsi  t  C 


»H "  Z'\r>Z'F'  riZ'F' 


0F»j  column:  0 


uj  column:  *mkt(CcCZ  Z'  +  VC?”  z'  )Ctf  *"  Cg 


AV'Z'^Z'F'  nZ'F' 


*i  column:  OttP.C.T*'  +  VC”''Z')R%  -  M)C,*,pC,f 


^krt„„(cccZz'  +  vcym  *  HCtfCVC, 

zr r<ZP r>P'r' \nP' F’  °KV 


W  Z'  \/rtZP /-iZFrtP' F1 


rzv c>zf (iv v  \nr- 
°pl  °ito  °r»»  )kF' 


°  H 


B,  column:  emfcpCfn,(£cC11^'  +  0C£ z  )(C^Cp7C: 

r,ZPrZFrP'Pl\TjP,P' 

°pf  Wo  W»i  P'F'o 


(8.3.4~23f) 


r> 

— 1+ 

duf 

+ 


119 


P/4  column:  0 


<j>ir  column:  0 
<j>y j  column:  0 
column:  0 
ii pit j  colunm:  0 
6pnj  column:  0 
iij  column:  0 

u'0  column:  -  eklm(CeC^'p'  +  VC%'P'  )CP'F' 


03  column:  -  tkii{CcCYk  P  +  VC ^  p  )<?£  F 


Pj4  column:  0 


<j>fr  column:  0 
<f>32  column:  0 
(j)^3  column:  0 


upnj  column:  0 
column:  0 


row: 


uj  column: 
u'p  column: 

03  column: 

P4  column: 

column: 
<^j2  column: 
$43  column: 


0 

[(CcC";'z'  +  VC** 

0 

0 

0 

0 

0 


)RQZ2  -  M)Cf? 


92k3 

du'^du'p 


u  p" j  column:  0 
dpiij  column:  0 


uj  column:  0 
u'p  column:  0 


(8.3.4-23h) 


(8.3.4— 23i) 


03  column:  0 


The  elements  of  E  and  F  are  determined  from  perturbations  of  Eqs.  (8.3.4-15)  which 
govern  the  lift,  drag,  and  pitching  moment.  Thus  the  E  matrix  may  be  defined  as: 


Cc  row: 


&zi  column:  paCciW§t  -f  \p«c~^Wr§i+ 


1  iWzi^q 

rp-c  i4G™ 


Wz2  column:  p.ccjH’^  -  jP„c^W'£1+ 


-7rpac 


wrG 


Q 

Z 12 


G§12  column:  ^pac2\W\ 


(8.3.4-24a) 


V  row: 

column:  paccdW§i  +  \PaC~^^Z2 
Wrz2  column:  pnccdW§2  ~  \p*c~^^z\ 


(8.3.4-24b) 


<?^12  column:  0 

Cnc  row: 

Wz\  column:  0 

.  0  (8.3.4-24c) 

W%2  column:  0 

G§j2  column:  0 
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M  row: 


W§ j  column:  p„c2cmW^,  +  \pac7~^^Zi' 

1  i  fi'S  >,0 

16  W‘  MG"12 

column:  pac2cmw£2  - 

1  c>"jLg« 

IS  Pa  \W\  212 


(8.3>4-24d) 


(?f12  column: - —  7r/ottc3  ITT^I 


and  for  the  elements  of  jP 


£r  row: 


Hr^j  column:  0 


:  Q 

W Z2  column:  0 


(8.3.4-25a) 


V  row: 


Gz  12  column:  0 


:  <? 

Tl^j  column:  0 
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W 22  column:  0 


(8.3.4-25b) 


Cnc  row: 
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Gz  12  column:  0 
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i  Q 

W Z2  column:  0 


(8.3.4-25c) 
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Gz  12  column:  — — 
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M  row: 


W zi  column:  — 

16 


:  Q  , 

W zi  column:  0 


(8.3.4~25d) 


,  3  npac* 

GZIJ  column:  - 

From  the  relations  defining  the  relative  velocity  components  and  gradient,  the  elements  of 
G  and  H  can  be  determined.  The  elements  of  G  are  defined  as 


WZa  row: 
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column:  0 
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u'0  column :  -  «*{m <?<f P  [ci\  A{U\  +  f ■ y A )  +  <?£  F  Vf,^]  — ~ 4 


s0^tmcr;F'  cf/nzing 


dKr 


du'i 


^Xc^c^rYac^RV  +  cfi*  Jjft 


dKr 


d*i 
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row: 
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<j)fr  column: 
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(j)A2  column: 

0 

4>a3  column: 

0 

up"j  column: 

+  c, 

(8.3.4~26b) 
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9p»j  column:  €jki  j  C3k  F  fi£,  ( 4- 

r Z'Aj>QF'  \/rZ'ArF'A  ,  rZ'ArF'A\_ 

— Cu  Kf, ,  +C/23  t/A3  j- 
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«Jj  column:  eumC33  Clk  ^F'*^r+ 
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The  elements  of  H  are  defined  as: 
W%a  row: 

Pj  column:  —  C^A 

•  A 

$lr  column:  0 


<^12  column:  —  C%lAR<%2A 
in  column:  -  C^ARQAA 


upi'j  column: 


dp»j  column: 


Uj  column: 
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(8.3,4-27a) 
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Gf12  row: 


-.A 

Pj  column:  0 


^lr  column:  0 


<j>i2  column:  -  A 

^13  column:  -  C^AC^A 
uf"  j  column:  0 


0F"j  column:  C{'j 


7,‘F' 


(8.3.4-27b) 


Uj  column:  0 


u'p  column: 


0k3 


63  column:  C33 

8.3.5.  Spatial  Discretization 

The  variables  u;  and  83  are  expanded  in  a  set  of  polynomials  based  on  reference  30. 
The  “CO”  functions  (u3  and  83)  are  expanded  in  terms  of  ^i(x)  where  *  =  *3 /£.  The 
functions  used  beyond  the  first  two  standard  linear  functions  are  orthonorinalized.  The 
Cl  functions  (ua)  are  expanded  in  terms  of  /?,•(*).  The  functions  used  beyond  the  first 
four  standard  cubic  functions  are  orthonormalized.  The  details  of  the  orthonormalization 
procedure  are  specified  below. 

The  expansions  are  given  by 

nq 

««  =  XI  9a»ft(*) 

t=] 

n9 

u3  =  »&(*) 

i=  1 


(8.3. 5-1) 


N* 

1=1 
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The  functions  0,  for  i  >  2  and  0i  for  i  >  4  are  constructed  from  the  Jacobi  polynomials 
gn(x)  =  Gn-i(p1q‘,x)  where  p  5  and  9  =  3  for  the  CO  functions  and  where  p  =  9  and 
9  =  5  for  the  Cl  functions.  # 

Letting  x  =  the  CO  shape  functions  are 


01  =1  —  S 


02  =* 


(8. 3. 5-2) 


0<  =®(1  -  x)gi-2(x)fi-2 

where  3  <  i  <  N  +  1  and  JV  =  N3  ot  N4.  The  recursion  relations  used  to  compute  the 
polynomials  are 
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0l(*)  =1 


where 


g2(X)  =X  -  - 

gj{x)  =g2{x)gj-i(x)  -  gj-2(x)Aj.2 

0i(*)  =° 
g'(x)  =1 

Sj(x)  =9j-iix)  +  92{x)g'j.1(x)  -  g'j_3(x)Aj-2 

g'l(x)  =0 
*"(»>  =o 

g'!{x)  =2<^_,(x)  +  g2{x)g'j.l{x)  - 

jlir]i(;  +  ,-i)(.+p-iK.-;;,rv).  p  =  M,3 

(2i  +  p  ~  2)(2t  +  p  -  1  )2 (2£  +  p) 

/,2  =30 


(8. 3.5-3) 


(8.3.5-4) 


^*2 =s 

The  derivatives  of  the  shape  functions  are  then 
*1  =  -  1 
<  =0 
i>\ ;  =i 
V>2  =0 

=((1  -  2x)«7i_2  +  *(1  "  ®)sJ_21/.-3 


(8. 3.5-5) 


V»i'  =I-25i-j(*)  -  4{x  -  -)^_2(*)  +  *(1  - 
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Similarly,  the  Cl  shape  functions  are 

/?,  =1  -  3x2  +  2x3 

02  =x  -  2x2  +  x3 

03  =3x2  -  2x3  (8.3. 5-6) 

04  =x3  —  x2 

Pj  =*2(1  -  ®)20«_4 (*)/«- 4 

where  5  <  t  <  N  +  1  and  N  =  jVj  or  N2  •  The  g' s  are  the  same  as  above  for  the  CO  shape 
functions,  .4,  is  evaluated  for  p  =  9,  q  =  5,  and  / 2  =  630.  The  higher  derivatives  are 

g'l'ix)  =g'2"(x)  =  0 

g'j'(x)  =3g'j_i(x)  +  ff2(z)5y_,  -  g'jL2{x)Aj.2 

(8.3.5-7) 

0l»  «*?"(*)  =  0 

»T(*)  +  ~  WMAj- 2 

and  derivatives  of  the  shape  functions  are 

0\  -r  -6x(l  -  x)  0';  ^  12(x  -  |)  /?"'  =  12  /?;»'  =  0 

$  =  3x2  -  4x  +  1  0%  =  6(x  -  §)  02  =  6  /?2 "  =  0  o  e_o\ 

$  =  6x(l  -  x)  /?"  =  - 12(x  -  i)  /?"'  =  -12  /?""  =  0 

0'4  -  3x2  -  2x  /?"  =  G(x-J)  01"  =  6 


#  - |-4x(l  -  x)(x  -  |)(/i_4(x)  4-  x2(l  -  x)2#j_4 (*)]/»•  .4 

4'  —  { 1 2(«2  -  x  I  | )^,_.*(x)  -  8x(l  -  x)(x  -  §)p-_4  +  x2(l  -  x)2flf-'i4(x)]/,_4 

0"'  -|24(x  -  i)gt-4{x)  +  3G(x2  -  x  +  |)fl'_4(x)  -  12x(l  -  x)(x  -  ±)g"_4{x)+ 
x2(l  -  x)2r/J'l4(x)]/,-4 


0"n  ^(24</, -4(j)  +  9G(x  -  |)</,'_4(a;)  +  ~2(x2  -  x  +  |K_4(x)- 
lGx(l  -  x)(x  -  !)(/!" 4(x)  +  x2(l  -  x)2(/J"'4(x)J/i_4 


(8. 3. 5-9) 


These  formulas  fur  shape  functions,  when  substituted  into  expressions  for  virtual  work 
of  either  internal,  inertial,  ur  applied  loads,  produce  integrands  that  depend  only  on  x3. 
These  integrals  tan  be  evaluated  to  any  accuracy  desired  by  use  of  Gaussian  quadrature. 
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8.3.6.  Transformation  from  Nodal  Coordinates 

In  GRASP,  a  different  set  of  generalized  coordinates  are  used  for  the  beam  element 
than  those  for  the  nodes.  It  is  therefore  necessary  to  calculate  the  beam  generalized 
coordinates  in  terms  of  the  nodal  displacements  and  rotational  variables  at  both  the  root 
and  tip  of  the  beam,  so  that  the  beam  equations  can  be  written  using  a  convenient  set  of 
generalized  coordinates. 

The  beam  generalized  coordinates  q„i  for  i  =  1,2, 3, 4  determine  the  uQ  displacements 
at  the  beam  root  and  tip.  Similarly,  gja  determines  the  U3  displacement  at  the  root  and 
tip,  and  q\a  determines  the  #3  rotation  at  the  root  and  tip.  The  exact  relations  are 

9.1  =ufl«' 


912 

914 

922 

924 


-Wi 


_ rT'R 

-°31 


_r.R,R 

_032 


~ °32 


(qn\ 


923 

932 


=cRr 


Uf 


(8. 3.6-1) 


<741  =  sin 


r 


942  —  sin 


cr'R  -.qT'TqTr 


U' 


The  rotation  expressions  are  derived  from  expressions  for  CP  "F ’  written  in  terms  of 
and  63  (see  ref.  24,  equations  4,  17,  and  60-62)  for  which 


Cy |  F  —  sin  62  —  Wj 


_  nil  r Ol  ,  .  | 

C32  —  —  cos  62  sinffi  =  u2 

CPi  F  =  —  cos  62  sin O3  —  —  yj  1  -  u', 2  sin  63 


(S.3.6-2) 
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8.3.7.  Transformation  io  Forces  and  Moments 

The  generalized  forces  calculated  for  the  beam  element  root  and  tip  correspond  to 
the  beam  generalized  coordinates.  These  forces  must  now  be  transformed  into  forces  and 
moments  at  the  root  and  tip  nodes.  The  virtual  work  at  the  root  can  be  written  in  terms  of 
the  static  residuals  Qr  and  the  linear  coefficient  matrices.  In  terms  of  the  beam  generalized 
coordinates,  this  relation  is 


SWn  =  6  qnr(-Qn  +  Lnqn) 


(8.3.7-1) 


where  6qnr  -  |6?ji  &/21  bq3]  Sqi2  bq22  6<7-nJ  and  =  [hi  hi  hi  O12  hi  94 1 J »  and  Lr 
is  a  linear  operator  representing  [Mnjjj  *  Cltjt  +  Rr)-  Note  that  this  equation  defines 
the  negative  of  the  virtual  work.  The  explanation  for  treating  the  virtual  work  in  this 
manner  is  that  it  is  conventional  for  Lr  to  he  positive,  and  Lr  is  normally  considered  to 
he  positive  on  the  left  hand  side  of  the  equations  of  motion,  while  Qn  is  positive  on  the 
right-lmnd  side. 

The  root  node  virtual  displacements  and  rotations  maj  he  related  tu  the  beam  virtual 
generalized  coordinates  by  the  expression 


'fyn ' 

Hq2\ 
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1  = 
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3xl 

H'n 

bq22 

3x3 
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3xl 

1  4 

V  ^94 1  - 

(8.3.7  2) 


where  the  root  node  virtual  displacements  are  6ur  -  6R ^  n ,  and  the  root  nude  virtual 
rotations  are  -  ^V'/j  11  ■  The  6x6  coefficient  matrix  that  prcmultiplics  the  root  node 
virtual  displacement  and  rotation  vector  .s  called  Tr  and  matrix  Rr  (ref.  24,  cq.  67)  is 


Rr 


0  -teg* 

teg*  0  teg* 

,.«*«  ^n'n 

()  -~2i_ —  -  -Las., — 

1  1 


(8. 3. 7-3) 


Similarly .  the  perturbed  mot  nude  displacements  u r  and  rotations  Or  arc  related  to 
the  perturbed  element  generalized  coordinates  qn  through  the  expression 


On 


(8.3.7  4) 
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When  the  virtual  work  at  the  beam  root  is  transformed  into  nodal  coordinates  by 
the  substitution  of  equations  (8.3. 7-2)  and  (8. 3. 7-4)  into  equation  (8.3.7-1),  the  following 
expression  is  obtained: 

~6Wn  =|£urT  6i/>rt}Trt  +  LRTR  j 

(8.3. 7-5) 

f>i>RT  J  |{-3rtq«}  +  trtlrtr 


First  consider  the  transformation  of  LR,  which  contains  the  dynamic  matrices  MRy 
C[ j,  and  Kr.  The  transformation  of  the  element  generalized  coordinates  into  the  nodal 
generalized  coordinates  introduces  the  transformation  matrix  TR  into  the  expression  for  the 
virtual  work.  Since  TR  is  a  function  of  CR  R,  which  is  a  function  of  the  nodal  rotations, 
it  must  also  be  perturbed  to  recover  any  additional  perturbation  contributions.  In  the 
case  of  the  linearized  dynamic  matrices  MRy  CR,  and  KR,  no  new  perturbation  terms  are 
introduced  by  the  transformation,  since  any  such  contributions  would  be  nonlinear. 

The  transformation  of  the  static  generalized  force  QR  is,  however,  another  matter.  In 
this  case,  transformation  does  contribute  an  additional  term,  called  the  geometric  stiffness 
term  K £,  to  the  linearized  dynamic  equations.  Geometric  stiffness  originates  from  the 
perturbation  of  TR. 

-GqrQr  =  ~[6urT  HrT\^'rT{Qr}  (8.3.7-6) 

When  TR  is  perturbed, 


where 


T  ci  Ti^Tr1 


-6<!rQr  =  -  6t/>rj  J 


=  -  L^U«T  J 


dgR 

0Trt 

dqR 


<irQr 


Tr 


{::) 


(8. 3. 7-7) 


Qr 
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dqR 
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3X3 
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(8. 3. 7-8) 


When  equation  (8.3.7  7)  is  multiplied  out,  only  one  of  the  3x3  submatrices  is  nonzero. 
This  submatrix  is  called  the  root  geometric  stiffness  matrix  k(jR,  and  it  contains  only  terms 
that  are  related  to  the  nodal  rotations. 


xg4  =  - 


BRrt 

dqn 


RrOrQro 


(8.3. 7-9) 
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where 

(8.3.7-10) 

Written  in 

index  notation  to  allow  the  isolation  of  koR,  the  root  geometric  stiffness  is 

KcRi; 0Rj  =  RR{i8R.QRek 

(8.3.7-11) 

» 

and 

kGn{i  =  d(jR  RRii  Qroh 

(8.3.7-12) 

The  geometric  stiffness  matrix  used  to  transform  all  of  the  root  nodal  degrees  of  freedom 
is  then 


0  0 
3x3  3x3 

0  kaR 

3X3  3x3 


(8.3.7-13) 


The  virtual  work  at  the  root  can  now  be  written  in  the  form 

-SWr  =[£uhT  8ipRT J  |-TnTC?H  +  TrTLrTr 

where 

Qr  =™rTQr 

L*r  =r.TRTLRTR  +  TRTK% 


(8.3.7-14) 


(8.3.7-15) 


The  transformation  of  the  generalized  forces  and  moments  at  the  tip  of  the  element 
into  nodal  forces  and  moments  is  similar  to  that  for  the  root.  In  beam  element  generalized 
coordinates,  the  virtual  work  at  the  tip  is 

— t  =  8qri'T  (— Qt  +  Lrqr)  (8.3.7-16) 

where  Sqr1  =  [8qp  8q23  8q32  8qu  8q24  j  and  q%  =  [qu  q23  q32  ql4  q24  q42 J,  and 

Lt  is  a  linear  operator  representing  +  CV~  +  KT ).  Note  the  similarity  with 

equation  (8. 3. 7-1). 
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The  equation  that  relates  the  tip  node  virtual  displacements  to  the  element  virtual 
generalized  coordinates  is 

'Sqi a' 

^923 

,  (m  , 

oqu 

6q2i 

>  ^942  - 

where  the  tip  node  virtual  displacements  are  and  the  tip  node  virtual 

rotations  are  Sip r  =  BiPrR-  The  6x6  coefficient  matrix  that  premultiplies  the  tip  node 
virtual  displacement  and  rotation  vector  is  called  Tt  and  matrix  Rt  is 
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(8.3.7-17) 


0 

ecf,* 

-icTS 
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l-(6’f'«p 

1  -ic%*y 

(8.3.7-18) 


Similarly,  the  perturbed  tip  node  displacements  ut  and  rotations  are  related  to 
the  perturbed  element  generalized  coordinates  qr  through  the  expression 


9r 


(8.3.7-19) 


The  expression  for  the  virtual  work  at  the  tip  is  similar  to  the  expression  for  the 
virtual  work  at  the  root. 

-SWT  =  [6utT  6il>rr J  {{~TtTQt}  +  TtTLtTt  (8.3.7-20) 


As  in  the  derivation  of  the  transformation  of  Lr ,  no  additional  terms  result  from  the 
transformation  of  L r •  There  art  however,  geometric  stiffness  terms  that  result  from  the 
transformation  of  Qt •  Following  the  derivation  of  the  root  geometric  stiffness, 


where 


jP 

-SqrQr  —  ~[Sutt  6iprT J  — T-  Ty  Qt  (8.3.7-21) 


dqr 


8Tt 

dqr 


0  0 

3x3  3x3 
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When  equation  (8.3.7-21)  is  expanded,  only  one,  nonzero  3x3  submatrix  remains.  It  is 
called  the  tip  geometric  stiffness  matrix  kGT>  and 


JcgtOt  =  -Ctr?~—RtCtr0tQt9 
dqT 

(8.3.7-23) 

where 

(8.3.7-24) 

Written 

in  index  notation  to  allow  the  isolation  of  &gt>  tl.e  tip  geometric  stiffness  is 

ho  t„  =  -Cl«a-^Rrm,C«jQT„, 

(8.3.7-25) 

or 

^GT{j  =  {CTRkGTCRT)ij 

(8.3.7-26) 

where 

kora  =  ^  RtgQto, 

(8.3.7-27) 

Therefore, 
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The  virtual  work  at  the  tip  can  now  be  written  in  the  form 
— Wt  =\5utt  6ij)TT J  j—  TttQt  +  TtTLtTt  |  4-  || 

=l««rT  «*rTj(- Q J  +  ir  {]£}} 

where 

Qt  =TttQt 


L*r  =TttLtTt  +  TttK$ 


(8.3.7-29) 


(8.3.7-30) 


For  both  the  root  and  tip,  derivatives  of  Rn  and  Rr  with  respect  to  the  qtJ  are  needed. 
The  only  nonzero  elements  of  these  arrays  may  determined  from 
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where 


where  C  i 
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CR'n  at  the  root  and  CT'n  at  the  tip  and  R  is  Rr  at  the  root  and  Tty  at  the 
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9.  CONCLUDING  REMARKS 


In  response  to  the  limitations  of  previous  methods  for  analyzing  rotorcraft,  GRASP 
has  been  developed,  GRASP  is  a  general-purpose  program  which  treats  the  nonlinear  static 
and  linearized  dynamic  behavior  of  rotorcraft  represented  by  arbitrarily  connected  rigid- 
body  and  beam  elements.  Large  relative  motions  and  deformation-induced  displacements 
and  rotations  are  permitted  (as  long  as  the  strains  in  the  beam  element  are  small).  Periodic 
coefficients  arc  not  treated,  restricting  the  solutions  to  rotorcraft  in  axial  flight  and  on  the 
ground. 

GRASP  uses  a  modern  approach  for  modeling  structures,  incorporating  the  features 
of  several  traditional  methods.  The  basic  approach  which  provides  the  foundation  for 
large  relative  motion  kinematics  is  derived  from  “multibody”  research  with  an  expanded 
emphasis  on  multiple  levels  of  substructures.  This  is  combined  with  the  finite  element 
approach  which  provides  flexible  modeling  through  the  use  of  libraries  of  elements,  con¬ 
straints,  and  nodes.  The  use  of  a  variable-order  polynomial  beam  element  makes  the  finite 
element  approach  more  effective.  The  incorporation  of  aeroelastic  effects,  including  inflow 
dynamics  and  nonlinear  aerodynamic  coefficients  for  the  beam  element,  further  extends 
the  capabilities. 

Due  to  the  fact  that  GRASP  was  developed  using  structured,  modular,  software  meth¬ 
ods,  changes  to  the  code  are  relatively  easy  to  perform.  This  makes  it  practical  to  modify 
the  code  in  order  to  enhance  its  functionality.  Some  of  the  many  areas  where  possibilities 
for  enhancements  exist  are  expanded  solution  procedures,  improved  aerodynamic  models, 
expanded  modeling  capabilities,  new  elements,  and  new  constraints. 

Existing  solution  procedures  (steady-state  and  asymmetric  eigenproblem)  could  easily 
be  expanded  to  include  a  symmetric  eigensolution.  This  solution  procedure  would  take  the 
symmetric  part  of  the  linearized,  perturbed  equations  of  motion,  then  calculate  the  eigen¬ 
values  and  eigenvectors.  The  symmetric  eigensoiution  would  be  to  generate  the  modes  for 
another  new  solution  procedure,  the  subspace  reduction.  The  subspace  reduction  would 
allow  the  user  to  solve  for  the  asymmetric  eigensolution  using  a  reduced  set  of  admissible 
functions.  A  reference  deformations  solution  procedure  would  allow  a  user  to  take  any 
steady  state  solution  and  use  it  either  as  an  initial  guess  for  another  steady-state  problem, 
or  ns  the  state  about  which  the  linearization  is  performed  for  an  eigensolution.  The  refer¬ 
ence  deformations  solution  would  lift  the  restriction  that  the  same  model  must  be  used  in 
the  the  steady  state  solution  and  the  eigensolution.  Another  valuable  enhancement  would 
be  to  extend  GRASP  to  forward  flight  using  either  a  time-domain  solution,  a  periodic 
solution,  or  both. 

Enhancements  to  the  aerodynamics  could  include  adding  the  capability  for  table- 
lookup  for  the  aerodynamic  coefficients,  and  perhaps  making  those  coefficients  functions 
of  Mach  number.  Another  possibility  would  be  to  incorporate  a  lifting  line  or  lifting  surface 
theory  to  calculate  the  aerodynamic  forces.  Wake  geometry  could  also  be  included.  Other 
valuable  enhancements  to  the  aerodynamic  model  would  be  the  inclusion  of  transonic  and 
dynamic  stall  effects. 
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The  modeling  capabilities  could  also  be  improved  with  the  addition  of  the  ability  to 
model  applied  loads.  It.  might  be  advantageous  to  include  simple,  dead  loads  (forces), 
and  geometrically  nonlinear  loads  such  as  applied  moments  and  nonconservative  forces. 
With  the  rapid  growth  of  control  theory,  some  sort  of  control  representation  should  be 
included  in  GRASP.  This  could  be  as  simple  as  specifying  the  thrust  level  of  the  rotor,  or 
as  complex  as  a  complete  control  representation  including  sensors,  actuators,  and  control 
laws.  In  addition,  it  would  be  convenient  to  implement  a  “generic”  node.  Such  a  node 
would  be  used  to  allow  the  user  to  define  generalized  coordinates  not  associated  with  any 
of  the  predefined  nodes. 

GRASP  would  greatly  benefit  from  the  addition  of  a  composite  beam  element  and 
a  direct-input  element.  The  composite  beam  element  would  be  able  to  rigorously  treat 
the  structural  couplings  introduced  by  composite  layups.  This  element  might  also  include 
the  effects  of  shear  deformation,  initial  curvature,  and  warping  rigidity.  The  direct-input 
element  would  be  used  in  conjunction  with  the  generic  node  to  allow  the  user  to  define  the 
properties  of  elements  that  are  not  included  in  GRASP.  An  example  of  such  a  use  would  be 
taking  a  set  of  inodes  from  a  NASTRAN  analysis  to  represent  the  fuselage  of  a  helicopter. 

New  constraints  that  would  enhance  the  capabilities  of  GRASP  include  a  moving- 
frame  constraint,  a  pin  constraint,  and  a  clamp  constraint.  The  moving-frame  constraint 
would  allow  a  frame  to  deform  with  the  structure.  Currently,  frame  motion  is  independent 
of  the  structure.  The  pin  constraint  would  allow  a  node  to  rotate  arbitrarily  about  either 
a  frame  or  another  node.  Eliminating  all  motion  of  a  node  would  be  accomplished  using 
the  clamp  constraint. 

From  this  description  of  possible  enhancements,  it  should  be  obvious  that.  GRASP 
has  a  great  potential  for  growth.  Because  of  its  modular  construction,  GRASP  has  the 
capability  to  handle  expansion  without  requiring  massive  rewriting  of  the  existing  equations 
and  code.  This  framework  makes  GRASP  a  desirable  platform  for  future  development. 
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